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Thermodynamics of type II superconductors in electromagnetic field based on the Ginzburg - 
Landau theory is presented. The Abrikosov flux lattice solution is derived using an expansion 
in a parameter characterizing the "distance" to the superconductor - normal phase transition 
line. The expansion allows a systematic improvement of the solution. The phase diagram of the 
vortex matter in magnetic field is determined in detail. In the presence of significant thermal 
fluctuations on the mesoscopic scale (for example in high T c materials) the vortex crystal melts 
into a vortex liquid. A quantitative theory of thermal fluctuations using the lowest Landau level 
approximation is given. It allows to determine the melting line and discontinuities at melt, as 
well as important characteristics of the vortex liquid state. In the presence of quenched disorder 
(pinning) the vortex matter acquires certain "glassy" properties. The irreversibility line and static 
properties of the vortex glass state are studied using the "replica" method. Most of the analytical 
methods are introduced and presented in some detail. Various quantitative and qualitative features 
are compared to experiments in type II superconductors, although the use of a rather universal 
Ginzburg - Landau theory is not restricted to superconductivity and can be applied with certain 
adjustments to other physical systems, for example rotating Bose - Einstein condensate. 
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I. INTRODUCTION 

Phenomenon of superconductivity was initially defined 
by two basic properties of classic superconductors (which 
belong to type I, see below): zero resistivity and perfect 
diamagnetism (or Meissncr effect). The phenomenon was 
explained by the Bose - Einstein condensation (BEC) 
of pairs of electrons (Cooper pairs carrying a charge 
— e* = — 2e, constant e* considered positive throughout) 
below a critical temperature T c . The transition to the su- 
perconducting state is described phenomenologically by a 
complex order parameter field ^ (r) = |^ (r)| e lx ^ with 
| 'J' | 2 proportional to the density of Cooper pairs and its 
phase x describing the BEC coherence. Magnetic and 
transport properties of another group of materials, the 
type II superconductors, are more complex. An external 
magnetic field H and even, under certain circumstances, 
electric field do penetrate into a type II superconductor. 
The study of type II superconductor group is importance 
both for fundamental science and applications. 

A. Type II superconductors in magnetic field 

1. Abrikosov vortices and some other basic concepts 

Below a certain field, the first critical field iJ cl , the 
type II superconductor is still a perfect diamagnct, but 
in fields just above H cl magnetic flux does penetrate the 



material. It is concentrated in well separated "vortices" 
of size A, the magnetic penetration depth, carrying one 
unit of flux 



The superconductivity is destroyed in the core of a 
smaller width £ called the coherence length. The type 
II superconductivity refers to materials in which the ra- 
tio n = A/£ is larger than k c = 1/^/2 (Abrikosov, 1957). 
The vortices strongly interact with each other, forming 
highly correlated stable configurations like the vortex lat- 
tice, they can vibrate and move. The vortex systems in 
such materials became an object of experimental and the- 
oretical study early on. 

Discovery of high T c materials focused attention to cer- 
tain particular situations and novel phenomena within 
the vortex matter physics. They are "strongly" type II 
superconductors n <~ 100 >> k c and are "strongly fluctu- 
ating" due to high T c and large anisotropy in a sense that 
thermal fluctuations of the vortex degrees of freedom are 
not negligibly, as was the case in "old" superconductors. 
In strongly type II superconductors the lower critical field 
H c i and the higher critical field H C 2 at which the mate- 
rial becomes "normal" are well separated H c2 /H cl ~ k 2 
leading to a typical situation H c i « H < H c2 in which 
magnetic fields associated with vortices overlap, the su- 
perposition becoming nearly homogeneous, while the or- 
der parameter characterizing superconductivity is still 
highly inhomogeneous. The vortex degrees of freedom 
dominate in many cases the thermodynamic and trans- 
port properties of the superconductors. 

Thermal fluctuations significantly modify the proper- 
ties of the vortex lattices and might even lead to its melt- 
ing. A new state, the vortex liquid is formed. It has 
distinct physical properties from both the lattice and the 
"normal" metal. In addition to interactions and ther- 
mal fluctuations, disorder (pinning) is always present, 
which may also distort the solid into a viscous, glassy 
state, so the physical situation becomes quite compli- 
cated leading to rich phase diagram and dynamics in 
multiple time scales. A theoretical description of such 
systems is a subject of the present review. Two ranges of 
fields, H « H c2 and H » H c \, allow different simplifi- 
cations and consequently different theoretical approaches 
to describe them. For large k there is a large overlap of 
their applicability regions. 



2. Two major approximations: the London and the 
homogeneous field Ginzburg - Landau models 

In the fields range H « H c2 vortex cores are well 
separated and one can employ a picture of line-like vor- 
tices interacting magnetically. In this approach one ig- 
nores the detailed core structure. The value of the order 
parameter is assumed to be a constant \T/o with an ex- 
ception of thin lines with phase winding around the lines. 
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FIG. 1 Schematic magnetic phase diagram of a type II super- 
conductor. 



Magnetic field is inhomogeneous and obeys a linearized 
London equation. This model was developed for low T c 
superconductors and subsequently elaborated to describe 
the high T c materials as well. It was comprehensively 
described in numerous reviews and books (Blatter et al., 
1994; Brandt, 1995; Kopnin, 2001; Tinkham, 1996) and 
will not be covered here. 

The approach however becomes invalid as fields of or- 
der of H C 2 are approached, since then the cores cannot be 
considered as linelike and profile of the depressed order 
parameter becomes important. The temperature depen- 
dence of the critical lines is sketched in Fig. 1. The 
region in which the London model is inapplicable in- 
cludes typical situations in high T c materials as well as 
in novel "conventional" superconductors. However pre- 
cisely under these circumstances different simplifications 
are possible. This is a subject of the present review. 
When distance between vortices is smaller than A (at 
fields of several H C 2) the magnetic field becomes homo- 
geneous due to overlaps between vortices. This means 
that magnetic field can be described by a number rather 
then by a field. This is the most important assumption of 
the Landau level theory of the vortex matter. One there- 
fore can focus solely on the order parameter field \& (r) . 
In addition, in various physical situation the order pa- 
rameter 1 F is greatly depressed compared to its maximal 
value ^o, due to various "pair breaking" effects like tem- 
perature, magnetic and electric fields, disorder etc. For 
example, in an extreme case of H ~ H C 2 only small "is- 
lands" between core centers remain superconducting, yet 
superconductivity dominates electromagnetic properties 
of the material. One therefore can rely on expansion 
of energy in powers of the order parameter, a method 
known as the Ginzburg - Landau (GL) approach, which 
is briefly introduced next. 

To conclude, while in the London approximation one 
assumes constant order parameter and operates with de- 
grees of freedom describing the vortex lines, in the GL 
approach the magnetic field is constant and one operates 
with key notions like Landau wave functions describing 



the order parameter. 



B. Ginzburg - Landau model and its generalizations 

An important feature of the present treatise is that 
we discuss a great variety of complex phenomena using 
a single well defined model. The mathematical methods 
used are also quite similar in various parts of the review 
and almost invariably range from perturbation theory to 
the so called variational gaussian approximation and its 
improvements. This consistency often allows to consider 
a smooth limit of a more general theory to a particular 
case. For example a static phenomenon is obtained as a 
small velocity limit of the dynamical one, the clean case 
is a limit of zero disorder and the mean field is a limit 
of small mesoscopic thermal fluctuations. The model is 
motivated and defined below, while methods of solution 
will be subject of the following sections. The complexity 
increases gradually. 



1. Landau theory near T c for a system undergoing a second 
order phase transition 

Near a transition in which the U (1) phase symmetry, 
vj/ — > e iX \l/, is spontaneously broken a system is effectively 
described by the following Ginzburg-Landau free energy 
(Mazenko, 2006): 



F[9] 



drdz{ 
a'\m 2 - 



2m* 
b' 



2m* 



|V*| 



(2) 



*| 4 } + Fn 



Here r = (x, y) and we assumed equal effective masses in 
the x — y plane m* = = m*, both possibly different 
from the one in the z direction m*/m* = 7^. This antic- 
ipates application to layered superconductors for which 
the anisotropy parameter j a can be very large. The 
last term, F n , the "normal" free energy, is independent 
on order parameter, but might depend on temperature. 
The GL approach is generally an effective mesoscopic ap- 
proach, in which one assumes that microscopic degrees of 
freedom are " integrates out" . It is effective when higher 
powers of order parameter and gradients, neglected in 
eq.(2) are indeed negligible. Typically, but not always, it 
happens near a second order phase transition. 

All the terms in eq.(2) are of order (1 — i) 2 , where 
t = T/T c , while one neglects (as "irrelevant") terms 
of order (1 — t) like \^f\ 6 and quadratic terms contain- 
ing higher derivatives. Generally parameters of the GL 
model eq.(2) are functions of temperature, which can be 
determined by a microscopic theory or considered phc- 
nomenologically. They take into account thermal fluctu- 
ations of the microscopic degrees of freedom ("integrated 
out" in the mesoscopic description). Consistently one 
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expands the coefficients "near", with coefficient a' van- 
ishing at T c as (1 — t): 



a'(T) 



a(l-t)+a' (l-tf + ... 



(3) 



b'(T) = b' + b" (l-t) + 
to* (T) = to* + to*' (1 — t) + . 



The second and higher terms in each expansion are omit- 
ted, since their contributions are also of order (1 — t) 3 
or higher. Therefore, when temperature deviates signifi- 
cantly from T c , one cannot expect the model to provide 
a good precision. Minimization of the free energy, eq.(2), 
with respect to * below the transition temperature deter- 
mines the value of the order parameter in a homogeneous 
superconducting state: 



1*1 



l*o| 2 (l-0, 



l*o| 



Substituting this into the last two terms in the square 
bracket in eq.(2), one estimates them to be of order 
(1-t) , while one of the terms dropped, I*! 6 , is indeed 
of higher order. The energy of this state is lower than 
the energy of normal state with \& — 0, namely, F n by 



^ = -F eL (l-t),F GL ^|*o| 4 



(5) 



is the condensation energy density of the superconductor 
at zero temperature. 

The gradient term determines the scale over which fluc- 
tuations are typically extended in space. Such a length 
£, called in the present context the coherence length, is 
determined by comparing the first two terms in the free 
energy: 

V 2 * ~ r 2 * oc (1 - t) % £ = j=^== . (6) 

1 /2 

So typically gradients are of order (1 — t) 1 , and the first 
term in the free energy, eq.(2) is therefore also of the or- 
der (1 — t) 2 . Since the order parameter field describing 
the Bose - Einstein condensate of Cooper pair is charged, 
minimal coupling principle generally provides an unam- 
biguous procedure to include effects of electromagnetic 
fields. 



2. Minimal coupling to magnetic field. 

Generalization to the case of magnetic field is a 
straightforward use of the local gauge invariance prin- 
ciple (or the minimal substitution) of electromagnctism. 
The free energy becomes 



F[*,A] 



I 



dvdz\ 



l 2to* 



|D*| 2 + /? 



urn 2 



b' 



1*1 



2to* 
G n [A] , 



(7) 



while the Gibbs energy is 



G [*, A] = F [*] + j 



(b-h)- 

8tt 



(8) 



Here B = V x A and we will assume that "external" 
magnetic field (considered homogeneous, see above) is 
oriented along the positive z axis, H =(0,0, H). The co- 
variant derivatives are defined by 



_ „ 2tt a 
D = V + i— A. 



(9) 



The "normal electrons" contribution G n [A] is a part 
of free energy independent of the order parameter, but 
can in principle depend on external parameters like tem- 
perature and fields. Minimization with respect to * and 
A leads to a set of static GL equations, the nonlinear 
Schrodinger equation, 

-£-G = — -— — D 2 ^ - ^—D 2 z t> + a'* + 6'|*| 2 * = 0, 

Sip* 2m* 2to* z 11 

(10) 

and the supercurrent equation,. 

S _ c „ „ 
c—G = — V x B — J s — J 
SA Air 



0, 



where the supercurrent and the "normal" current 

ie*h 



Js = 



2to 

ie*h 
2^* 



_ (**D* - *D**) 



(11) 



(12) 



(#*V* - W**) - —A |*| 2 

cm* 



J„ can be typically represented by the Ohmic conductiv- 
ity J„ = cr„E, and vanishes if the electric field is absent. 

Comparing the second derivative with respect to A 
term in eq.(ll) with the last term in the supercurrent 
equation eq.(12), one determines the scale of typical vari- 
ations of the magnetic field inside superconductor, the 
magnetic penetration depth: 



V 2 A~A- 2 (l-t)A 
This leads to 



=*2 



,lg-A(l-i)|*o| 2 . (18) 



mi 



2e* V tt< i r, 



(14) 



The two scales' ratio defines the GL parameter k = A/£. 
The second equation shows that supercurrent in turn is 
small since it is proportional to \^\ 2 < *o. Therefore 
magnetization is much smaller than the field, since it 
is proportional both to the supercurrent creating it and 
to 1/k 2 . Since magnetization is so small, especially in 
strongly type II superconductors, inside superconductor 
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B w H and consistently disregard the " supercurrent" 
equation eq.(ll). Therefore the following vector poten- 
tial 



A = {-By,0,0)~(-Hy,0,0) 



(15) 



(Landau gauge) will be use throughout. The validity of 
this significant simplification can be then checked apos- 
teriori. 

The upper critical field will be related in section II to 
the coherence length eq.(6) by 



Hc2 — 



_$o_ 

2tt£ 2 



(16) 



The energy density difference between the superconduc- 
tor and the normal states Fgl in eq.(2) can therefore be 
reexpressed as 



F ( 



GL 



H c2 



(17) 



3. Thermal fluctuations 



Thermal fluctuations on the microscopic scale have al- 
ready been taken into account by the temperature de- 
pendence of the coefficients of the GL free energy. How- 
ever in high T c superconductors temperature can be high 
enough, so that one cannot neglect additional thermal 
fluctuations which occur on the mcsoscopic scale. These 
fluctuations can be described by a statistical sum: 



Z = Jv$> (r) X>** (r) exp j- 



F [**,*] 



(18) 



where a functional integral is taken over all the config- 
urations of order parameter. In principle thermal fluc- 
tuations of magnetic field should be also considered, but 
it turns out that they arc unimportant even in high T c 
materials (Dasgupta and Halpcrin, 1981; Halperin et al, 
1974; Herbut and Tesanovic, 1996; Herbut, 2007; Lobb, 
1987) . 

Ginzburg parameter, the square of the ratio of T c to the 
superconductor energy density times correlation volume, 



Gi = 2 



l^FGLiHc 



& 



(19) 



generally characterizes the strength of the thermal fluc- 
tuations on the mesoscopic scale (Ginzburg, 1960; Larkin 
and Varlamov, 2005; Levanyuk, 1959) and where $o = 
^ . The definition of Gi is the standard one as in (Blatter 
et al, 1994) contrary to the previous definition used early 
in our papers, for example in (Li and Rosenstein, 2002a, 
2003). Here £ c = 7 a " 1 £ is the coherence length in the 
field direction. The Ginzburg parameter is significantly 
larger in high T c superconductors compared to the low 
temperature one. While for metals this dimensionless 
number is very small (of order 10 -6 or smaller), it be- 
comes significant for relatively isotropic high T c cuprates 



like YBCO (10~ 4 ) and even large for very anisotropic 
cuprate BSCCO (up to Gi = 0.1 — 0.5). Physical reasons 
behind the enhancement arc the small coherence length, 
high T c and, in the case of BSCCO, large anisotropy 
7 a <~ 150. Therefore the thermal fluctuations play a 
much larger role in these new materials. In the presence 
of magnetic field the importance of fluctuations is further 
enhanced. Strong magnetic field effectively suppresses 
long wavelength fluctuations in direction perpendicular 
to the field reducing dimensionality of the fluctuations 
by two. Under these circumstances fluctuations influence 
various physical properties and even lead to new observ- 
able qualitative phenomena like the vortex lattice melt- 
ing into a vortex liquid far below the mean field phase 
transition line. 

Several remarkable experiments determined that the 
vortex lattice melting in high T c superconductors is first 
order with magnetization jumps (Bcidcnkopf et al, 2005, 
2007; Nishizaki et al, 2000; Willemin et al, 1998; Zeldov 
et al, 1995), and spikes in specific heat (it was found that 
in addition to the spike there is also a jump in specific 
heat which was measured as well) (Bouquet et al, 2001; 
Lortz et al, 2006, 2007; Nishizaki et al, 2000; Schilling 
et al, 1996, 1997). These and other measurements like 
the resistivity and shear modulus point towards a need to 
develop a quantitative theoretical description of thermal 
fluctuations in vortex matter (Liang et al, 1996; Matl 
et al, 2002; Pastoriza et al, 1994) To tackle the difficult 
problem of melting, the description of both the solid and 
the liquid phase should reach the precision level below 1% 
since the internal energy difference between the phases 
near the transition temperature is quite small. 



4. Quenched Disorder. 

In any superconductor there are impurities either 
present naturally or systematically produced using the 
proton or electron irradiation. The inhomogeneities both 
on the microscopic and the mesoscopic scale greatly af- 
fect thermodynamic and especially dynamic properties 
of type II superconductors in magnetic field. Abrikosov 
vortices are pinned by disorder. As a result of pinning 
the flux flow may be stopped and the material restores 
the property of zero resistivity (at least at zero tempera- 
ture, otherwise thermal fluctuations might depin the vor- 
tices) and make various quantities like magnetization be- 
comes irreversible. Disorder on the mesoscopic scale can 
be modeled in the framework of the Ginzburg - Landau 
approach adding a random component to its coefficients 
(Larkin, 1970). The random component of the coefficient 
of the quadratic term W(r) is called 5T disorder, since 
it can be interpreted as a local deviation of the critical 
temperature from T c . The simplest such a model is the 
"white noise" with local variance: 



a' -> a' [1 + W (r)] ; W (r) W (r') = n^&S (r - r') . 



(20) 
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A dimensionless disorder strength n, normalized to the 
coherence volume, is proportional to the density of the 
short range point - like pinning centers and average 
"strength" of the center. The disorder average of a static 
physical quantity A, denoted by "A" in this case, is a 
gaussian measure p [W] 

A = J VW{r)A[W]p[W], (21) 

p[W) = Ne ,7V- 1 = / VW (r) e wHc . 

The averaging process and its limitations is the subject 
of section IV, where the replica formalism is introduced 
and used to describe the transition to the glassy (pinned) 
states of the vortex matter. They are characterized by 
irreversibility of various processes. The quenched disor- 
der greatly affects dynamics. Disordered vortex matter 
is depinned at certain "critical current" J c and the flux 
flow ensues. Close to J c the flow proceeds slowly via 
propagation of defects (elastic flow) before becoming a 
fast plastic flow at larger currents. The I-V curves of the 
disordered vortex matter therefore are nonlinear. Dis- 
order creates a variety of "glassy" properties involving 
slow relaxation, memory effects etc. Thermal fluctua- 
tions in turn also greatly influence phenomena caused by 
disorder both in statics and dynamics. The basic effect 
is the thermal depinning of single vortices or domains of 
the vortex matter. The interrelations between the inter- 
actions, disorder and thermal fluctuations are however 
very complex. The same thermal fluctuations can soften 
the vortex lattice and actually can also cause better pin- 
ning near peak effect region . Critical current might have 
a "peak" near the vortex lattice melting. 

C. Complexity of the vortex matter physics 

In the previous subsection we have already encoun- 
tered several major complications pertinent to the vortex 
physics: interactions, dynamics, thermal fluctuations and 
disorder. This leads to a multitude of various "phases" 
or states of the vortex matter. It resembles the com- 
plexity of (atomic) condensed matter, but, as we will 
learn along the way, there are some profound differences. 
For example there is no transition between liquid and 
gas and therefore no critical point. A typical magnetic 
(T — B) phase diagram advocated here(Li et al, 2006b) 
is shown on Fig. 2b. It resembles for example, an exper- 
imental phase diagram of high T c superconductor (Di- 
vakar et al, 2004; Sasagawa et al, 2000) LaSCO Fig. 
2a. Here we just mention various phases and transi- 
tions between them and direct the reader to the rele- 
vant section in which the theory can be found. Let us 
start the tour from the low T and B corner of the phase 
diagram in which, as discussed above, vortices form a 
stable Abrikosov lattice. Vortex solid might have several 
crystalline structures very much like an ordinary atomic 
solid. In the particular case shown at lower fields the 




FIG. 2 Magnetic phase diagram of high T c . (a) Experimen- 
tally determined phase diagram of LaSCO(Divakar et al., 
2004). (b) Theoretical phase diagram advocated in this arti- 
cle. 

lattice is rhombic, while at elevated fields in undergoes a 
structural transformation into a square lattice (red line 
on Fig. 2). These transitions are briefly discussed in 
section II. Thermal fluctuations can melt the lattice into 
a liquid (the "melting" segment of the black line), sec- 
tion III, while disorder can turn both a crystal and a 
homogeneous liquid into a "glassy" state, Bragg glass or 
vortex glass respectively (section IV). The correspond- 
ing continuous transition line (blue line on Fig. 2) is often 
called an irreversibility line since glassiness strongly af- 
fects transport properties leading to irreversibility and 
memory effects. 

To summarize we have several transition lines 

1) The first order (Bouquet et al, 2001; Schilling et al, 
1996, 1997; Zeldov et al, 1995) melting line due to ther- 
mal fluctuations was shown to merge with the "second 
magnetization peak" line due to pinning forming the 
universal order - disorder phase transition line (Fuchs 
et al, 1998; Radzyner et al, 2002). At the low tempera- 
tures the location of this line strongly depends on disor- 
der and generally exhibits a positive slope (termed also 
the "inverse" melting (Paltiel et al, 2000a, b)), while in 
the " melting" section it is dominated by thermal fluctua- 
tions and has a large negative slope. The resulting max- 
imum at which the magnetization and the entropy jump 
vanish is a Kauzmann point (Li and Rosenstein, 2003). 
This universal "order - disorder" transition line (ODT), 
which appeared first in the strongly layered supercon- 
ductors (BSCCO (Fuchs et al, 1998)) was extended to 
the moderately anisotropic superconductors (LaSCCO 
(Radzyner et al, 2002)) and to the more isotropic ones 
like YBCO (Li and Rosenstein, 2003; Pal et al, 2001, 
2002). The symmetry characterization of the transition 
is clear: spontaneous breaking of both the continuous 
translation and the rotation symmetries down to a dis- 
crete symmetry group of the lattice. 

2) The "irreversibility line" or the "glass" transition 
(GT) line, which is a continuous transition (Deligiannis 
et al, 2000; Senatore et al, 2008; Taylor et al, 2003; 
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Taylor and Maple, 2007). The almost vertical in the T—B 
plane glass line clearly represents effects of disorder al- 
though the thermal fluctuations affect the location of 
the transition due to thermal depinning. Experiments 
in BSCCO (Beidenkopf et al, 2005, 2007; Fuchs et al, 
1998) indicate that the line crosses the ODT line right 
at its maximum, continues deep into the ordered (Bragg) 
phase. This proximity of the glass line to the Kauzmann 
point is reasonable since both signal the region of close 
competition of the disorder and the thermal fluctuations 
effects. In more isotropic materials the data are more 
confusing. In LaSCCO (Divakar et al, 2004; Sasagawa 
et al, 2000) the GT line is closer to the "melting" section 
of the ODT line still crossing it. It is more difficult to 
characterize the nature of the GT transition as a "sym- 
metry breaking" . The common wisdom is that " replica" 
symmetry is broken in the glass (either via "steps" or 
via "hierarchical" continuous process) as in the most of 
the spin glasses theories (Dotsenko, 2001; Fischer and 
Hertz, 1991). The dynamics in this phase exhibits zero 
resistivity (neglecting exponentially small creep) and var- 
ious irreversible features due to multitude of metastable 
states. Critical current at which the vortex matter starts 
moving is nonzero. It is different in the crystalline and 
homogeneous pinned phases. 

3) Sometimes there are one or more structural tran- 
sitions in the lattice phase (Divakar et al, 2004; Es- 
kildsen et al, 2001; Gilardi et al, 2002; Jaiswal-Nagar 
et al, 2006; Johnson et al, 1999; Keimer et al, 1994; 
Li et al, 2006a; McK. Paul et al, 1998; Sasagawa et al, 
2000). They might be cither first or second order and 
also lead to a peak in the critical current (Chang et al, 
1998a, b; Klironomos and Dorsey, 2003; Park and Husc, 
1998; Rosenstein and Knigavko, 1999). 



D. Guide for a reader. 

1. Notations and units 

Throughout the article we use two different systems 
of units. In sections not dealing with thermal fluctua- 
tions, namely in section II and section IVA we use units 
which do not depend on "external" parameters T and H, 
just on material parameters and universal constants (for 
example the unit of length is the coherence length £). 
More complicated parts of the review involving thermal 
fluctuations utilize units dependent on T and H. For 
example the unit of length in directions perpendicular to 

the field direction becomes magnetic length I = i\J~~^ff-- 
However throughout the review basic equations and im- 
portant results, which might be used for comparison with 
experiments and other theories, are also stated in regular 
physical units. 



a. The mean field units and definitions of dimensionless pa- 
rameters Ginzburg - Landau free energy, eq.(2), con- 



tains three material parameters m* , m* (in the a—b direc- 
tions perpendicular to the field and in the field direction 
respectively), aT c ,b'. If in addition the ST C disorder, 
introduced in cq.(20), is present, it is described by the 
disorder strength n. These material parameters are usu- 
ally expressed via physically more accessible lengths and 
time units £, £ c , A. 



h 



^/2m* c aT c 



(22) 



Despite the fact that one often uses temperature depen- 
dent coherence length and penetration depth, which as 
seen in equation eqs.(6) and (13) might be considered as 
divergent near T c , we prefer to write factors of (1 — t) 
explicitly. 

From the above scales can form the following dimen- 
sionless material parameters Gi and 



K = V£ 5 7a = m* c /m* 



(23) 



From the scales one can form units of magnetic and 
electric fields, current density and conductivity: 



H, 



<:2 



$0 

2< 



2 ' 



(24) 



as well as energy density Fql- These can be used to de- 
fine dimensionless parameters, temperature T, magnetic 
and electric fields H, E 



t = 



b = 



B 

H C 2 



h = 



H 

H C 2 ' 



(25) 



from which other convenient dimensionless quantity de- 
scribing the proximity to the mean field transition line 
arc formed 



a H 



l-t-b 



(26) 



The unit of the order parameter field (or square root of 
the Cooper pairs density) is determined by the mean field 
value |* | 2 = 



* = 



V2|* | 



b' 



2aT r 



1/2 



(27) 



and the Boltzmann factor and the disorder correlation 
in the physics units (length is in unit of £ in x — y plane 
and in unit of £ c along c axis, order parameter in unit as 
defined by the equation above) is 



F 



T 



-I / d?x{\\DM 2 + 
Ut J £ 



\\d^\ 2 

G[*,A] 
T 



1 - t 



{l + W(r))\i/>\ 



W (r) W (r') 



F [**,*] 
T 

nS (r 



— I d 3 x 



1 

r') ,u t = 



(b-h) 2 

4 

s/2Gint. 
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b. The LLL scaled units When dealing with thermal fluc- 
tuations, the following units depend on parameters T, 
H and E. Unit of length in directions perpendicular to 
the field can be conveniently chosen to be the magnetic 
length, 



(28) 



in the field direction, while in the field direction it is 
different: 



Gitb 



-1/3 



(29) 



Motivation for these fractional powers of both tempera- 
ture and magnetic field will become clear in section III. 
we rescale the order parameter to ip by an additional 
factor: 



* = * 



Gitb 



1/3 







1>. 



(30) 



Instead of cih or o,h,e it will be useful to use "Thouless 
LLL scaled temperature" : (Ruggeri and Thouless, 1976; 
Ruggcri, 1978; Thouless, 1975) 



l-t-b 

/ , — s 2/3 ' 
/ VGitb ] ' 



(31) 



4 J 



The scaled energy is defined by 



'Gitb 



4/3 



2ttk 2 

and magnetization by 



/(or) , 



(32) 



M 

Hc2 



1 



Gitb 



2/3 



Atth 2 



m (ax) ■ 



(33) 



The disorder is characterized by the ration of the strength 
of pinning to that of thermal fluctuations 



_ (i-*r 

r vGiVH n - 



(34) 



techniques. We chose to describe some of them in de- 
tail, while others are just mentioned in the last section. 
We do not describe numerous results obtained using the 
elasticity theory or numerical methods like Monte Carlo 
and molecular dynamics simulations, although compari- 
son with both is made, when possible. 

The techniques and special topics include: 

1) Translation symmetries in gauge theories (electro 
- magnetic translations) in IIA. Their representations, 
the quasi - momentum basis (IIIB) is used throughout 
to discuss excitations of vortex matter either thermal or 
elastic. 

2) Perturbation theory around a bifurcation point of a 
nonlinear PDE (differential equations containing partial 
derivatives). This is very different from the perturbation 
theory used in linear systems, for example in quantum 
mechanics 

3) Variational gaussian approximation to field theory 
(Klcinert, 1995) is widely used in III to IV. It is de- 
fined in IIIC in the path integral form and subsequently 
shown to be the leading order of a convergent series of 
approximants, the so called optimized perturbation se- 
ries (OPE). The next to leading order, the post gaus- 
sian approximation, is related to the Cornwall -Jackiw 
-Tomboulis method is sometimes used, while higher ap- 
proximants are difficult to calculate and are obtained to 
date for the vortex liquid only. 

4) Ordinary perturbation theory in field theory is de- 
veloped in the beginning of every section with enough 
details to follow. Spatial attention is paid to infrared 
(IR) and sometimes ultraviolet (UV) divergencies. We 
generally do not use the renormalization group (RG) re- 
summation, except in subsection HID, where it is pre- 
sented in a form of Borcl - Pade approximants. 

5) Replica method to treat quenched disorder is intro- 
duced in IVB and used to describe the static and the ther- 
modynamic properties of pinned vortex matter. Most of 
the presentation is devoted to the replica symmetric case, 
while more general hierarchial matrices are introduced in 
IVD following Parisi's approach (Mezard, 1991; Parisi, 
1980). 

Some technical details are contained in Appendix. We 
compare with available experiments on type II supercon- 
ductors in magnetic field, while application or adaptation 
of the results to other fields in which the model can be 
useful (mentioned in summary) are not attempted.. 



3. Results 



2. Analytical methods described in this article 

Discussion of properties of the GL model in magnetic 
fields utilizes a number of general and special theoretical 



All the important results (in both regular physical 
units and the special units described above) are provided 
in a form of Mathematica file, which can be found on our 
web site. 
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II. MEAN FIELD THEORY OF THE ABRIKOSOV 
LATTICE 

In this section we construct, following Abrikosov orig- 
inal ideas (Abrikosov, 1957), a vortex lattice solution of 
the static GL equations eq.(10) "near" the H C 2 (T) line. 
In a region of the magnetic phase diagram in which the 
order parameter is significantly reduced from its maxi- 
mal value eq.(4), one does not really see well sepa- 
rated "vortices" since, as explained in the previous sec- 
tion, their magnetic fields strongly overlap. Very close to 
H C 2 (T) even cores approach each other and consequently 
the order parameter is greatly reduced. Only small "is- 
lands" between the core centers remain superconducting. 
Despite this superconductivity dominates electromag- 
netic, transport and sometimes thermodynamic proper- 
ties of the material. One still has a well defined "centers" 
of cores: zeroes of the order parameter. They still repel 
each other and thereby organize themselves into an or- 
dered periodic lattice. 

To see this we first employ a heuristic Abrikosov's ar- 
gument, based on linearization of the GL equations and 
then develop a systematic perturbative scheme with a 
small parameter - the " distance" from the H C 2 (T) line 
on the T — H plane. The heuristic argument naturally 
leads to the lowest Landau level (LLL) approximation, 
widely used later to describe various properties of the 
vortex matter. The systematic expansion allows to as- 
certain how close one should stay from the H C 2 line in 
order to use the LLL approximation. Having established 
the lattice solution, spectrum of excitations around it 
(the flux waves or phonon) are obtained in the next sub- 
section. This in turn determines elastic, thermal and 
transport properties of vortex matter. 



A. Solution of the static GL equations. Heuristic solution 
near H c2 

1. Symmetries, units and expansion in kT 2 

Broken and unbroken symmetries 

Generally, before developing (sometimes quite elabo- 
rate) mathematical tools to analyze a complicated model 
described by free energy eq.(2) and its generalizations, 
it is important to make full use of various symmetries 
of the problem. The free energy (including the external 
magnetic field) is invariant under both the three dimen- 
sional translations and rotations in the x—y (a—b) plane. 
However some of the symmetries in the x — y plane are 
broken spontaneously below the H C 2(T) line. The sym- 
metry which remains unbroken is the continuous trans- 
lation along the magnetic field direction z. As a result 
the configuration of the order parameter is homogeneous 
in the z direction ^> (r, z) = ^> (r), r = (x, y). Hence the 
gradient term can be disregarded and the problem be- 
comes two dimensional (here we consider the mean field 
equations only, when thermal fluctuations or point - like 
disorder are present the simplification is no longer valid.). 



Units, free energy and GL equations 

To describe the physics near H C 2{T), it is reasonable 
to use the coherence length £ = h/ \J2m*aT c as a unit 



of length (assuming for simplicity to* 



to*) and 



the value of the field V&o at which the " potential" part is 
minimal, eq.(4), (times ^/2) will be used as a scale of the 
order parameter field 



x = x/£, y = y/£; * 



( b ' 



\2aT c 



1/2 



(35) 



while the (zero temperature energy) density difference 
between the normal and the superconductor states Fql 
of eq.(17) determines a unit of energy density. Therefore 
dimensionless 2D energy F = 8L ^f G l ' wnere L z is the 
sample's size in the field direction, and cq.(8) takes a 
form: 



F 



dxdy 



(b-h) 



(36) 

Dimensionless temperature and magnetic fields are t = 
T/T c , b = B/H c2 , h = H/H c2 and k = A/£. The units 
of temperature and magnetic field are therefore T c and 

H C 2 



$0 



The linear operator H is defined as 



H = -UDl + dl + b). 



(37) 



It coincides with the quantum mechanical operator of a 
charged particle in a constant magnetic field. The covari- 
ant derivative (with all the bars omitted from now on) is 
D x = d x — iby and the constant is defined as 



a H = 



l-t-b 



(38) 



It is positive in the superconducting phase and vanishes 
on the H C 2 (T) line, as will be shown in the next subsec- 
tion. The reason why H is "shifted" by a constant —b/2 
compared to a standard Hamiltonian of a particle in mag- 
netic field will become clear there. In rescaled units the 
GL equation takes a form: 

ffW- a H H+ |*| 2 \i> = 0. (39) 
The equation for magnetic field takes a form 



iCEijdjb = Di$ + ex. 



(40) 



with boundary condition involving the external field h. 
Expansion in powers of kT 2 

In physically important cases one is encounters 
strongly type II superconductors for which k » 1. For 
example all the high T c cuprates have k of order 100, and 
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even low T c superconductors which are useful in applica- 
tions have k of order 10. In such cases it is reasonable to 
expand the second equation in powers of «~ 2 : 



b = h + n~ 2 b^ + ...; 
* = * (0) +«:- 2 * (1) + .... 



(41) 



It can be seen from eqs.(39) and (40) that to leading or- 
der in k~ 2 magnetic field b is equal to the external field 
h considered constant. Therefore one can ignore eq.(40) 
and use external field in the first equation. Corrections 
will be calculated consistently. For example magnetiza- 
tion will appear in the next to leading order. 

From now on we drop bars over W and consider the 
leading order in n~ 2 . Even this nonlinear differential 
equation is still quite complicated. It has an obvious 
normal metal solution \I> = 0, but might have also a 
nontrivial one. A simplistic way to find the nontrivial one 
is to linearize the equation. Indeed naively the nonlinear 
term contains the " small" fields ^ compared to one in the 
linear term. This assumption is problematic since, for 
example the coefficient of the ^ term is also small, but 
will follow this reasoning nevertheless leaving a rigorous 
justification to subsection B. 



2. Linearization of the GL equations near H C 2. 

Naively dropping the nonlinear term in eq.(39), one is 
left with the usual linear Schrocdingcr eigenvalue equa- 
tion of quantum mechanics for a charged particle in the 
homogeneous magnetic field 



(42) 



The Landau gauge that we use, defined in eq.(15), still 
maintains a manifest translation symmetry along the x 
direction, while the y translation invariance is "masked" 
by this choice of gauge. Therefore one can disentangle 
the variables: 



x f(y)- 



(43) 



resulting in the shifted harmonic oscillator equation: 



1 



f = 



1-t 



-f. 



(44) 



where Y = k x /b is the y coordinate of the center of the 
classical Larmor orbital. For a finite sample k x is dis- 
cretized in units of while the Larmor orbital center 
is confined inside the sample —L y /2 < Yt; < L y /2 lead- 
ing to BL ^ V = Nl values of k x . 

Nontrivial / (y) =/= solutions of the linearized equa- 
tion exist only for special values of magnetic field, since 
the operator H has a discrete spectrum 



for any Y (the Landau levels are therefore Nl times de- 
generate). These fields b^ satisfy 



1 



(46) 



and the eigenfunctions are: 



<f>Nk x (r) 



-1/4 



xe 



\k x x- 



2 N Nl 



N 



6 1/2 (V - k x /b)} (47) 



where Hn (x) are Hermit polynomials. As we will see 
shortly, the nonlinear GL equation eq.(39) acquires a 
nontrivial solution also at fields different from 6/v. The 
solution with N = (the lowest Landau level or LLL, 
corresponding to the highest b N — 1) appears at the bi- 
furcation point 



1 - t - b {t) = 



(48) 



or a H = 0. It defines the H c2 (T) = H c2 (1 - T/T c ) line. 

For yet higher fields the only solution of nonlinear GL 
equations is the trivial one: ^> = 0. This is seen as 
follows. The operator H is positive definite, as its spec- 
trum eq.(45) demonstrates. Therefore for an < all 
three terms in the free energy eq. (36) arc non - nega- 
tive and in this case the minimum is indeed achieved by 
^ = 0. For an > the minimum of the nonlinear equa- 
tions should not be very different from a solution of the 
linearized equation at B = H C 2 (T). 

Since the LLL, B = H c2 (T), solutions 



(49) 



are degenerate, it is reasonable to try the most general 
LLL function 



*W = E^^(r) 



(50) 



as an approximation for a solution of the nonlinear GL 
equation just below H c2 (T). However how should one 
chose the correct linear combination? Perhaps the one 
with the lowest nonlinear energy: the quartic term in en- 
ergy eq.(36) will lift the degeneracy Unfortunately the 
number of the variational parameters in eq.(50) is clearly 
unmanageable. To narrow possible choices of the coeffi- 
cients, one has to utilize all the symmetries of the lattice 
solution. Therefore we digress to discuss symmetries in 
the presence of magnetic field, the magnetic translations, 
returning later to the Abrikosov solution equipped with 
minimal group theoretical tools. 



3. Digression: translation symmetries in gauge theories 



E N = Nb 



(45) 



Translation symmetries in gauge theories 
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FIG. 3 Symmetry of the vortex lattice. Unit cell. 



Let us consider a solution of the GL equations invariant 
under two arbitrary translations vectors. Without loss of 
generality one of them dican be aligned with the x axis. 
Its length will be denoted by d. The second is determined 
by two parameters: 



dx = d(l,0); d 2 = d(p,p') 



(51) 



We consider only rhombic lattices (sufficient for most ap- 
plications), which are obtained for p = 1/2. The angle 9 
between diand d 2 is shown on Fig. 3. Flux quantization 
(assuming one unit of flux per unit cell) will be 



d 2 p'b = 2tt; p' — - tan 



(52) 



For the hexagonal (also called sometimes triangular) FLL 
C\ = iCo = iC. Geometry and the flux quantization 
gives us now £ 2 e? A = which becomes (in rescaled 

units of £) 



d7 



A 



2tt 2 



(56) 



We are therefore left again with just one variational pa- 
rameter 



V3d, 



1)2 ; 



(57) 



0(1(1 



4 E 



s/5d/ 



Naive nonmagnetic translation in the " diagonal" direc- 
tion, see Fig. 3, now gives 

x + —,y-\ - — \=ie d ^ tp A (x,y) (58) 

This is again a "regauging" , which generally accompanies 
a symmetry transformation. The "magnetic translation" 
now will be 



Generally an arbitrary translation in the x direction in 
the particular gauge that we have chosen, eq.(15), is very 
simple 

T dl *(x,y)=*(x + d,y) = e^ d * (x, y) , (53) 

where p = — iV is the "momentum" operator. Periodic- 
ity of the order parameter in the x direction with lattice 
constant d (in units of £, as usual) means that the wave 
vector k x in eq.(49) is quantized in units of =?■: k x = ^f-ro, 
n = 0, ±1, ±2, ...and the variational problem of eq.(50) 
simplifies considerably: 



*(r) = 



(54) 



^n(r) 



b 1 / 2 



Periodicity with lattice vector d2 is only possible only 
when absolute values of coefficients |C„| are the same 
and, in addition, their phases are periodic in n. 
Hexagonal lattice 

In this case the basic lattice vectors are di = d/\ (1,0), 
d 2 = d A (1/2, s/3/2) , see Fig. 3, 9 = 60°. As a next 
simplest guess to construct a lattice configuration out of 
Landau harmonics one can try a two parameter Ansatz 
Cn+2 = C n : 

V(x,y) = Co ^/'b-^e^-^-^Xhh) 

n even 

+C\ V 7r -l/4 6 -l/2 e ^™^-|(3/-|^i») 2 _ 
n odd 



T d2 = e 

The normalization is 
1 

vol 



A 2/ gH^P«i 2 Py) , (59) 



\<Po(x,y)\ =1- 



(60) 



which gives: \c\ 2 = S 1 / 4 ^ 1 / 2 /?. Combining the even and 
the odd parts, the normalized function also can be writ- 
ten in a form 



¥>A(r) = <^(& 1/2 r) 



(61) 



(r) = 3 1 / 8 J2 el[¥ 



2+3 1 /4 7r l/2 fa l_l(j / _ 3 l/4 7r l/2; ) 2 



l = — OO 



This form will be used extensively in the following sec- 
tions. 

General rhombic lattice 

All the rhombic lattices with magnetic field b are ob- 
tained from the Ansatz C n +2 = C n by assuming the 
phase C\ = iCo- 



tp(x,y,b) = J - — -= 2^ e u e 2 ^ . 

V ■ * I — DO 

(62) 

The hexagonal lattice corresponds to 9 = 60°, see Fig. 
3. One can check that a rhombic lattice indeed is invari- 
ant under magnetic translations by di and d 2 . The flux 
quantization takes a form 



1 ,2 „ 2tt 
-d e tan(9 = — . 



(63) 
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One notices d$ = dg(b) = de(l)/y/b,a,nd that generally 
we have a following relation, 



<p(x,y,b) = ip(b^ 2 x,b^ 2 y) 



(64) 



where the right side equation ip(x, y) is the solution in 
the case of b = 1 and we replace x, y by Vbx, Vby. There 
are of course infinitely many invariant functions differing 
by a " fractional" translation as well as by rotation of the 
lattice. These symmetries are "broken spontaneously" 
by the lattice. According to Goldstone theorem, they 
lead to existence of soft phonon modes in the crystalline 
phase and will be studied in section III. 

General magnetic translations and their alge- 
bra 

Let us generalize the discussion by considering an arbi- 
trary Landau gauge. Using the experience with regaug- 
ing of the two nontrivial translations in our gauge, which 
generally is defined a matrix 



B 






b 









(65) 



Magnetic translation operator for a general vector d 
should be defined as 



j> d — e -i(h d i B n+ r i B n) d j e id -P — e id '^, (66) 
with a generator P defined by 

Pi = -idi - ByiTj =Pi- Bjifj. (67) 
This can be derived using the general formula 

(68) 



e K £ L _ e K+L+\[K.L] 



valid when commutator [K, L] is proportional to the 
identity operator. Applying the formula to the case of 
the expression eq.(66) with K = —i (J^diBij + TiBij) dj, 
L = ip ■ d, and using the basic algebra [ri,Pj] = iSij, one 
indeed obtains a number 



[K, L] = [nBijdj^p ■ d] = idiBijdj. 



(69) 



The expression for magnetic translations can also be 
derived from a requirement that they commute with 
" Hamiltonian" H defined in eq.(37). In fact they com- 
mute with both covariant derivatives D i7 



Di — di + iBi-jrj, 



(70) 



as well. However, using the same basic algebra, one 
also observes that magnetic translations generally do not 
commute: Td x Td 2 differs from Td 2 T dl by a phase. This 
is a consequence of the Campbell-Bakcr-Hausdorff for- 
mula e K e L = e L e K e^ K ' L \ which follows immediately 
from cqs.(68) and (69) 



e »di-P e »d 2 -P _ e id 2 P e idi P e -[di P,d 2 P] 



(71) 



with the constant commutator given by [di • P, d2 • P] = 
ibdi x d 2 . The group property therefore is 



T dl T d2 = e - jbdlXd2 T d2 T dl , 



(72) 



from which the requirement to have an integer number 
of fluxons per unit cell of a lattice follows: 



6di x d 2 = 2ir x integer. 



(73) 



Note that the generator of magnetic translations is not 
proportional to covariant derivative Di = di — iB^rj. 
The relation is nonlocal, 



Pi iDi -\- ^ij^j, 



(74) 



where is the antisymmetric tensor. 



4. The Abrikosov lattice solution: choice of the lattice structure 
based on minimization of the quartic contribution to energy 

The Abrikosov (3 constant of a lattice structure 

To lift the degeneracy between all the possible "wave 
functions" with arbitrary normalization on the ground 
Landau level, one can try to minimize the quartic term in 
free energy eq.(36). It is reasonable to assume that more 
"symmetric" configurations will have an advantage. In 
particular lattices will be preferred over "chaotic" inho- 
mogeneous ones. Moreover hexagonal lattice should be 
perhaps the leading candidate due to its relative isotropy 
and high symmetry. This configuration is preferred the in 
London limit (Tinkham, 1996) since vortices repel each 
other and try to self assemble into the most homogeneous 
configuration. A simpler square lattice was considered in 
fact as the best candidate by Abrikosov and we start from 
this lattice to try to fix the variational parameter C. The 
energy constrained to the LLL is 



G 



I 



rfr 



-a H |vI/| 2 + I|*| 4 + ^(b-hr 



(75) 



The quartic contribution to energy density is propor- 
tional to the space average of \ip\ 4 which is called the 
Abrikosov /?<>: 



* - y 



do/ 2 f d o/2 

dx I dy\tp (x,y)\ 

d<>/2 



(76) 



= 12 dx d vJ2 



^-^{n\—ni+m—ni)x 



_b 

xe 2 



d<>/2 J-d^/2 
[(t/-doni) 2 +(j/-do™2) 2 + (j/-don3) 2 + (j/-do«4) 2 ] 



In principle one can slightly generalize the method we 
used to calculate analytically both integrals and sums 
(Saint- James et ai, 1969), however will refrain from do- 
ing so here, since in Appendix A a more efficient method 
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will be presented. The result is f3$ = 1.18. More gener- 
ally one it is shown there that for any lattice this constant 
is given by 



fcX"(»l,» 2 ) 



(77) 



is over the lattice 
For example for /?a = 



sites 
1.16 



where the summation 
X (711,712) =nidi + n 2 d 2 
for the hexagonal lattice. 

Energy, entropy and specific heat 

Free energy density of the leading order solution is in- 
deed negative. Substituting a variational solution, one 
has 



l-t 



b , ,2 

- M 



(78) 



+\\C\"\^\ i ] = -\C\"a H + \\C\ A ^. 



The FLL and the transition to the normal state can there- 
fore described well by a " dimensionally reduced" D = 
U(l) symmetric model with the complex "order parame- 
ter" C. It is similar to the Meissner state in the absence 
of magnetic field but in D = with the only difference 
being that between /3a and 1 (which is just about 10%). 
One minimizes it with respect to C: 



\C\ 2 = a H /0o. 



(79) 



The average energy density at minimum (still on the sub- 
space of square lattices) is given by 



— h 

vol J r 



2/3o 



(i-b-ty 



(80) 



or, returning to the unsealed units, the energy density is 



F 
vol 



H c2 



Attk 2 /3<> 



(81) 



The first derivative with respect to temperature T, the 
entropy density 



S= - 



H c2 



4itk 2 I3aT c 



(82) 



smoothly vanishes at transition to the normal phase. On 
the other hand the second derivative, the specific heat 
divided by temperature, jumps to a constant 



H c2 



8ttk 2 /3 a T 2 



(83) 



from zero in the normal phase. Note that in this section 
we use a simple GL model which neglects the normal 
state contribution to free energy, cq.(2), retaining only 
terms depending on the order parameter. The additional 
term is a smooth "background", also referred to as a 
" contribution of normal electrons" . 



Of course a similar argument is valid for any lattice 
symmetry with corresponding Abrikosov parameter (3 a- 
What is the correct shape of the vortex lattice? To min- 
imize the energy in this approximation is equivalent to 
the minimization of [3e with respect to shape of the lat- 
tice. This is achieved for the hexagonal lattice, although 
differences are not large. The square lattice incidentally 
has the largest energy among all the rhombic structures, 
some 2% higher than that of the hexagonal lattice. This 
sounds rather small, but for a comparison, the typical 
latent heat at melting (difference in internal energies be- 
tween lattice and homogeneous liquid) is of the same or- 
der of magnitude. 

Magnetization to leading order in 1/k 2 
Magnetization can be obtained via minimization of the 
Gibbs free energy with respect to magnetic induction B. 
In our units and within LLL approximation one can dif- 
ferentiate eq.(75) and the Maxwell term with respect to 
b: 



[h-b(r)] =47TK 2 m(r) = |*(r) 



(84) 



The magnetization m (r) is therefore proportional to the 
superfluid density |^ (r) | 2 and is thus highly inhomoge- 
neous. Its space average is 



m : 



J>(r) _C 2 \vo (r) 



Airvol 



Attk 2 



l-t-b 



(85) 



For large k (typical value for high T c superconductors is 
k = 100) the magnetization is of order 1/k 2 compared to 
H and therefore negligible. This justifies an assumption 
of constant magnetic induction, which can be slightly 
corrected: 



l-t 

'Wo 



K 2 h 



2/3 



l-t-h 



(86) 



Rescaling back to regular units, one has 



4tt v ' Auk 2 /?<> ' 



with 



1J T 

2 V Tr. 



B 

H C 2 



5" 



T H \ 

T c H C 2 J 



(87) 



, (88) 



valid up to corrections of order k~ 2 . 

A general relation between the current density 
and the superfluid density on LLL 

The pattern of supercurrent flow around vortex cores 
can be readily obtained by substituting the Abrikosov 
vortex approximation into the expression for the super- 
current density eq.(12). We derive here a general relation 
between an arbitrary static LLL function eq.(50) and the 
supercurrent. It will be helpful for understanding of the 
mechanism behind the flux flow, occurring in dynamical 
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situations, when electric field is able to penetrate a su- 
perconductor. The covariant derivatives acting on the 
LLL basis elements give: 



D, 



D 



= (d x - iby) {tt" i/%V*jK*-Uv-K/b? j 

= -i^^b 1 ' 2 (by - k x ) e^^lfe-^/ 6 ) 2 

= i{b/2) 1/2 4> N=1Mx (89) 

= -^- 1 /4 6 l/2 (6y _ fe) e ifex-|(y-fc,/6) 2 

= (b/2) 1/2 4>N=1,K. 



The covariant derivatives, which are linear combina- 
tions of "raising" and "lowering" Landau level operators, 



D x = iJ^ (S f +a) ; D y = J~^ (ot 



(90) 



-ic\, + d y - by 



-id x + d y + by 



2b 



2b 



therefore raise an LLL function to the first LL. One can 
check that the following relation is valid: 



i^f* (r) A* (r) + c.c. = e^dj 



*(r)|" 



(91) 



where is the antisymmetric tensor. We therefore have 
established an exact relation between the current density 



and (scaled with Jgl 
luid 

•Mr) 



c* 



to eq.(24)) superfluid density, 



Jgl 



■j according 



Ji (r) 



1 - 



Jgl 



a,(|*(r)| 2 ), (92) 



valid, however, on LLL states only. In regular units the 
current density is related to (unsealed) order parameter 
field by 



e*H 
2m*' 



^•(|*(r)| 2 ), (93) 



The supercurrent indeed creates a vortex around a dip 
in the superfluid density, Fig. 4. The overall current is of 
course zero, since the bulk integral is transformed into a 
surface one. An approximate solution described in this 
subsection is perhaps valid "near" the H c2 (T) line, how- 
ever to estimate the range of validity and to obtain a 
better approximation, one would prefer a systematic per- 
turbative scheme over an uncontrollable variational prin- 
ciple. This is provided by the an expansion. 





FIG. 4 Superflow around the vortex centers in the hexagonal 
lattice 



the development of the bifurcation point perturbation 
theory for the GL equation eq.(39). This type of the 
perturbation theory is quite different from the one used 
in linear equations like Schrodinger equation. 

One develops a perturbation theory in small an around 
the H c2 (T) line : 



* = Van h {0) + + a 2 H ¥ 2) 



(94) 



Note the fractional power of the expansion parameter 
in front of the "regular" series. This is related to the 
mean field critical exponent for a 4 type equation being 
1/2, so that all the terms in the free energy have the 
same power a 2 H and are "relevant", as we mentioned in 
Introduction. Substituting this series into eq.(39) one 

1/2 

observes that the leading (a^ ) order equation gives the 
lowest LLL restriction already motivated in the heuristic 
approach of the previous subsection: 



H^ (0) = 0, 



(95) 



resulting in VE^ ) = C^tp with normalization undeter- 
mined. It will be determined by the next order. The 
next to leading (a# 2 ) order equation is: 

- C<°V + C (0) |C (0) f V M 2 - 0. (96) 

Multiplying it with ip* and integrating over coordinates, 
one obtains 



C (0) 



<p\<p\ 



0. (97) 



B. Systematic expansion around the bifurcation point. 

1. Expansion and the leading order 

We have defined the operator H in eq.(37) in such a 
way that its spectrum will start from zero. This allows 



The first term vanishes since Hermitian operator H in 
the scalar product, defined as 



L x L y 



f* (r)g(r) = [f\g], 



(98) 
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can be applied on it and vanished by virtue of eq.(95). 
This way one recovers the "naive" result of eq.(79): 



C<°> 



1 



L X Ly J r 



C (0) 



/3 A = 0. (99) 



Note that to this order different lattices or in fact any 
LLL functions are " approximate solutions" . 



2. Higher orders corrections to the solution 
Next to leading order 

Higher order corrections would in principle contain 
higher Landau level eigenfunctions in the basis of solu- 
tions of the linearized GL equation eq.(42) for eigenval- 
ues En, eq.(47). As on the LLL for higher Landau levels 
one can combine them into a lattice with a certain (here 
hexagonal) symmetry: 



Pan (r) = / C ksc 4>Nk x (r) = <Pn ( fcl/2r ) ! 



(100) 



Vn (r) 



e «(^+3 1 / 4 7r 1 / 2 x)-|(8/-3 1 / 4 7r 1 / 



I=-c 



The coefficients are the same as given in the previous 
subsection, eq.(57). 

The order (a#-) l+1 ^ 2 correction can be expanded in the 
Landau levels basis, eq.(100) as 

oo 

¥ (r) = C®<p (V /2 r) + £ C$<p N (V /2 r) (101) 

JV=1 

(to simplify notations the LLL coefficient is denoted sim- 
ply C {i) rather than suppressing N = 0, the conven- 
tion we have been using already for tp = ipo). Inserting 

ux_:__ i. j 3/2. 

(102) 



this into eq.(96), one obtains to order a 3 ^ 2 

jr NbCtftpN = c (0 V - c (0) |c (0) | 

N=l 

The scalar product with tp N determines 



VM 2 



r w - 

<~- N - - 



Nbf3 3 / 2 ' 



where 



Pn = j \tp\ 2 (p N <P* 

vol /, 



(103) 



(104) 



To find C' 1 ' we need in addition also the order ajj 2 equa- 



tion: 



oo 

H»£ 2 = NbC N )( PN = *i-(C (0) ) 2 (2^ 1 |^| 2 + *^ 

JV=1 



^ 2 ). 

(105) 
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FIG. 5 Convergence of the bifurcation perturbation theory 



Inner product with ip gives 
3 



(106) 



The expansion can be relatively easily continued. Fig. 5 
shows three successive approximation. The convergence 
is quite fast even as far from the H c2 (T) line at for b = 
0.1, t = 0.5. 

Orders a 2 H and a A H in the expansion of free en- 
ergy. 

The mean field expression for the free energy to or- 
der a 2 H was obtained already using heuristic approach, 
eq.(80).. Inserting the next correction eqs.(106) and 
(103) into eq.(36) one obtains the free energy density: 



i?(2) + p(3) 

4 vol A 



'2(3 A 



Pib 



E 

JV=1 



(Pn? 
N 



(107) 



= -0.43a^- 0.0078^-, 



where A = „ c \ is the unit of energy density. 

It is interesting to note that Pn ^ only when n = 6j, 
where j is an integer. This is due to hexagonal symme- 
try of the vortex lattice (Lascher, 1965). For n — 6j it 
decreases very fast with j: 6 = -0.2787, /3 12 = 0.0249. 
Because of this the coefficient of the next to leading order 
is very small (additional factor of 6 in the denominator) . 
We might preliminarily conclude therefore that the per- 
turbation theory in an works much better that might be 
naively anticipated and can be used very far from tran- 
sition line. If we demand that the correction is smaller 
then the main contribution the corresponding line on the 
phase diagram will be b — 0.015 • (1 — t). For example 
the LLL melting line corresponds to an ~ 1 • This overly 
optimistic conclusion is however incorrect as calculation 
of the following term shows. 

How precise is LLL? 
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Now wc discuss in what region of the parameter space 
the expansion outlined above can be applied. First of 
all note that all the contributions to 'J'i are proportional 
to 1/6. This is a general feature: the actual expansion 
parameter is a#/6. One can ask whether the expansion is 
convergent and, if yes, what is its radius of convergence. 
Looking just at the leading correction and comparing it 
to the LLL one gets a very optimistic estimate. For this 
purpose higher orders coefficients were calculated (Li and 
Rosenstein, 1999a). The results for the \&2 are following: 



a 



(2) 
N 



1 

A6 



1 



(108) 



]T c$ (2[j\t,o|m,o] + [o,o|m,jv]); 



M=0 



and 



C (2 



3 Ft 



2/3a 



2^ 



(109) 



]T C£>C${[0,0\L,M\ + 2[M,0\L,0]), 



L,M=0 



where 



[K,L\M,N] 



1 f s|s * 



(110) 



We already can see that Cffi and are proportional 
to and in addition there is a factor of 1/N. Since, 



due to hexagonal lattice symmetry all the N ^= 6j 

(2) 

vanish, so do C N . We have checked that there is no 
more small parameters, so we conclude that the leading 
order coefficient is much larger than first (factor 6-5), 
but the second is only 6 times larger than the third. The 
correction to free energy density is 



4 vol A 



0.056 a\j 



(111) 



Accidental smallness by factor 1/6 of the coefficients in 
the an/b expansion due to symmetry means that the 
range of validity of this expansion is roughly a# < 66 
or B < Hc \^ ■ Moreover additional smallness of all the 
HLL corrections compared to the LLL means that they 
constitute just several percent of the correct result in- 
side the region of applicability. To illustrate this point 
we plot on Fig. 5 the perturbatively calculated solution 
for 6 = 0.1, t = 0.5. One can see that although the 
leading LLL function has very thick vortices (Fig. 5a), 
the first nonzero correction makes them of order of the 
coherence length (Fig. 5b). Following correction of the 
order (a/f/6) makes it practically indistinguishable from 
the numerical solution. Amazingly the order parameter 
between the vortices approaches its vacuum value. Para- 
doxically starting from the region close to H c2 the pertur- 
bation theory knows how to correct the order parameter 



so that it looks very similar to the London approxima- 
tion (valid only close to il cl ) result of well separated 
vortices. 

We conclude therefore that the expansion in au/b 
works in the mean field better that one can naively ex- 
pect. 



III. THERMAL FLUCTUATIONS AND MELTING OF 
THE VORTEX SOLID INTO A LIQUID 

In this section a theory of thermal fluctuations and 
of melting of the vortex lattice in type II superconduc- 
tors in the framework of Ginzburg - Landau approach is 
presented. Far from H C \(T) the lowest Landau level ap- 
proximation can be used. Within this approximation the 
model simplifies and results depend just on one param- 
eter: the LLL reduced temperature. To obtain an accu- 
rate description of both the vortex lattice and the vortex 
liquid different methods are applied. In the crystalline 
phase basic excitations are phonons. Their spectrum and 
interactions are rather unusual and the low temperature 
perturbation theory requires to develop a certain tech- 
nique. Generally perturbation theory to the two loop 
order is sufficient, but for certain purposes (like finding 
a spinodal in which metastable crystalline state becomes 
unstable) a self consistent "gaussian" approximation is 
required. In the liquid state both the perturbation the- 
ory and gaussian approximations are insufficient to get a 
precision required to describe the first order melting tran- 
sition and one utilizes more sophisticated methods. Al- 
ready gaussian approximation shows that the metastable 
liquid state persists (within LLL) till zero temperature. 
The high temperature renormalized series (around the 
gaussian variational state) supplemented by interpola- 
tion to a T = metastable "perfect liquid" state are 
sufficient. The melting line location is determined and 
magnetization and specific heat jumps along it are cal- 
culated. The magnetization of liquid is larger than that 
of solid by 1.8% irrespective of the melting temperature, 
while the specific heat jump is about 6% and decreases 
slowly with temperature. 



A. The LLL scaling and the quasi - momentum basis 
1. The LLL scaling 

Units and the LLL scaled temperature 

If the magnetic field is sufficiently high, we can keep 
only the N = LLL modes. This is achieved by enforcing 
the following constraint, 



ft 2 fte* „ T 
-D ^ = 



(112) 



2m* 2m*c 
where covariant derivatives were defined in eq.(9). Using 
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it the free energy eq.(8) simplifies: 

G [*, A] = / dr{^ \&M 2 +aT c (1 - t - b) |*| 
J 2m* 



(113) 



Originally the Ginzburg - Landau statistical sum, 
eq.(18), had five dimensionless parameters, three ma- 

1/2 

terial parameters k = A/£,7 a = (m*/m*) , and the 
Ginzburg number, defined by 



Gi 



e* 2 K 2 £T cla 
2irc 2 H 2 



(114) 



and two external parameters t = T/T c and b — B/H c2 . 
However, since there is now no gradient term in directions 
perpendicular to the field, it is missing one independent 
parameter. The Gibbs energy, 

e = -Tlogj^ex P -I Jg[*,B] J, (115) 

thus possesses the "LLL scaling" (Lee and Shenoy, 1972; 
Ruggeri and Thouless, 1976; Ruggcri, 1978; Thouless, 
1975). To exhibit these scaling relations, it is useful 
to use units of coordinates and fields, which are de- 
pendent not just on material parameters (as those used 
in section II), but also on external parameters, mag- 
netic field and temperature. One uses the magnetic 
length rather than coherence length as a unit of length 
in directions perpendicular to magnetic field, x = -y^x, 

y = -^y, while in the field direction different factor is 
used, z — ^- ^^J^) z. Magnetic field is rescaled 



as before with H c2 



while the order parameter field has 

2/3 

ip 2 . The use- 



/Gitb 



an additional factor: ^> 2 = 2^q 

fulness of the fractional powers additional factors will 
become clear later. 

The dimensionless Boltzmann factor becomes 



g[1>A 



G[*,A] 



i 



25/2^ 



4/3 



(116) 



+ 



-I 



k 2 ( VGitb \ r 
'^n ^ 4 J J 



2 5 / 2 
dr 



dr 



(b-hy 



l -\( k 4 > \ 2 + a T \M 2 + \\M i 



(117) 



where the LLL scale "temperature" is 



'Gitb 



a H = -- 



l-t-b 



'Gitb 



(118) 

The constant an was defined in eq.(38) and extensively 
used in the previous section. The scaled temperature 



therefore is the only remaining dimensionless parameter 
in eq.(116) in addition to the coefficient of the last term. 
Factors of 2 5 / 2 7r in definition of "dimensionless free en- 
ergy" / in cq.(116) are traditionally kept and will appear 
frequently in what follows. Assuming nonfluctuating con- 
stant magnetic field, one can disregard the last term in 
eq.(116), and consider the thermal fluctuations of the or- 
der parameter only. This assumption is typically valid 
in almost all applications and will be discussed in sub- 
section E. Certain physical quantities, the "LLL scaled" 
ones, are functions of this parameter only. We list the 
most important such quantities below. 

Scaled quantities 

The scaled free energy density is: 



fdM = — ^log j * 



(H9) 



where V is the rescaled volume and f(ar) is related to 
the free energy density in unsealed units by 

Turning to magnetization, let us return to conventional 
units eq.(113) and neglect fluctuations of magnetic field 
(considered in (Dasgupta and Halperin, 1981; Halpcrin 
et ai, 1974; Herbut and Tesanovic, 1996; Herbut, 2007; 
Lobb, 1987)). Within LLL magnetization in the presence 
of thermal fluctuations is determined from 



_6_ 

SB 



g = z~ 1 ( 4g[*,%^ /g| * !B1 -0. (121) 

J* SB 



Taking the derivative, one obtains 



-z- 1 



B-H 



47T 



(122) 



where from now on (...) denotes the thermal average. 
The magnetization on LLL is therefore proportional to 
the superfluid density 



M 



n C 2 



(123) 



This motivates the definition of the LLL scaled magneti- 
zation proportional to (l^l 2 ), 



m (ax) = — 



which is related to magnetization by 

2/3 



M 

H C 2 



1 



'Gi 



47TK 2 



tb 



m (or) . 



(124) 



(125) 
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Consequently 



M 



depends on ar only, the state- 



ment called "the LLL scaling" proposed in (Ruggeri and 
Thouless, 1976; Ruggeri, 1978; Tesanovic et al, 1992; 
Tesanovic and Andrccv, 1994; Thouless, 1975). It has 
been experimentally demonstrated in numerous experi- 
ments. 

The specific heat contribution due to the vortex matter 
is generally defined by C = -T-§^G{T,H). Usually, 
since the GL approach is applied near T c , one can replace 
T by T c in the Boltzmann factor, leaving the temperature 
dependence just inside the coefficient of \^\ 2 in eq.(113). 
In this case the normalized specific heat is defined as 



C 



a 



mf 



(126) 



where C m f = g^/t^lr 2 is the mean field specific heat 
of solid calculated in the previous section. Substituting 
G(T,H), if very near phase transition temperature, we 

can put t = 1 in the scaling factor ^-^-tb^n this case, we 
obtain: 



c= ~^h fd{aT) - 



(127) 



Since the range of applicability of LLL can extend be- 
yond vicinity of T c , especially at strong fields (since they 
depress order parameter), one should use a more compli- 
cated formula which docs not utilize T ~T r : 



16 — V3 Hb-l-t) {bt)2/3 



if 2 



xf(a T ) + 



(2-26 + t) 2 ,, 



9t 2 



(128) 



It no longer possesses the LLL scaling. 



2. Magnetic translations and the quasi - momentum basis 

It is necessary to use the representations of transla- 
tional symmetry in order to classify various excitations of 
both the Abrikosov lattice and a homogeneous state cre- 
ated when thermal fluctuations become strong enough. 
As we have seen in subsection IIB, presence of magnetic 
field makes the use of the translational symmetry a non- 
trivial task, due to the need to "regauge". Here we use 
an algebraic approach to construct the quasi - momenta 
basis and then to determine the excitation spectrum of 
the lattice and the liquid, which in turn determines its 
clastic and thermal properties. 

The quasi - momentum basis 

We motivated the definition of the magnetic transla- 
tion symmetries eq. (66) by the property that they trans- 
form various lattices onto themselves. More formally the 
x — y plane translation operators Ta, cq.(66), represent 
symmetries since they commute with " Hamiltonian" H of 
cq.(37). Excitations of the lattice are no longer invariant 



under the symmetry transformations. This in particular 
means that we cannot longer consider the problem as two 
dimensional. However, as in the solid state physics, it is 
convenient to expand them in the basis of eigenfunctions 
of the generators of the magnetic translations operators 
defined in eq.(66) and simple translations in the field, z, 
direction: 



~PfNk = k(p Nk ; T d ip Nk 
Pz<PNk = k z (p Nk ; 



e <PNk, 



(129) 



J-d z fNk — e ifNk- 



with commutation relation cq.(72): T^T^ = 
e~ tbdlXd2 Td 2 Td 1 ■ The tree dimensional quasi - momen- 
tum vector is denoted by k = (k, k z ). It is easy to con- 
struct these functions explicitly. On the N th Landau 
level the 2D quasi - momentum k function is given by: 



VNk (r) = Trip N (r) 



(130) 



where ki = Sijkj for i = x,y and (p^ (r) for a given lattice 
symmetry was constructed in IIA. Here we will take the 
hexagonal lattice functions of eq.(100). Indeed 



TdtpNk = T d T^tp N = e 1 x T~T A ip N = e l ' ifiNk- 

(131) 

To write it explicitly, the most convenient form of the 
magnetic translation is that of eq.(66), which gives 



(132) 



Since Td is unitary, the normalization is the same as that 
of (p j\[ . On LLL in our gauge one has: 

oo 

fk = e lxk *p (r + k) = 3 1 / 8 ( 133 ) 

/ — — CO 

e i [ 2lf +3 U^U2 (x+ky)l+xkx] _l [y _ kx _ 3 l/4 7T l/2 l] 2 

In the direction along the field one uses the usual mo- 
mentum: 



Vk (r) 



_ ik z z 



Vk (r) , 



(134) 



where, as before, we use the notation r = (r, z). 

The values of the quasi - momentum cover a Brillouin 
zone in the x — y plane. As usual, it is convenient to work 
in basis vectors of the reciprocal lattice, k = fcidi + fc 2 d 2 , 
with the basis vectors 



di = 



'3* 



V5F 

The measure is 



1, - 



1 

V5 



;d 2 = 



33 'o, 2 



v V3 



/ dk x dk y = 2tt dki / dk 2 ; 
Jb.z. Jo Jo 

f d 3 k = f dk z ( dk. 

J k J —oo J B.z. 



(135) 



(136) 
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Beyond LLL the quasi - momentum basis consists of 
ip^ (r), N th Landau level "wave functions" with quasi- 
momentum k: 



3V4 



2 N M ^ 

I — — CO 



H n {y-k x -&IWH) (137) 



xe 1 2 



ii + 3l/4 7r l/2 (x+feH) ; +xfeai] _l [y _ fe!c _3l/4 7r l/2; ] 2 



The construction is identical to LLL. Even in the homo- 
geneous liquid state, which is obviously more symmetric 
than the hexagonal lattice, we find it convenient to use 
this basis: 

1 r °° 

M = W^f2 / E (138) 

Energy in the quasi - momentum basis 

As was discussed in section II, the lowest energy con- 
figurations belong to LLL. There is an energy gap to 
any N > configuration, so it is reasonable that, for 
temperatures small enough, their contribution is small. 
Restricting the set of states over which we integrate to 
LLL 



xj){r) 



1 



(27T) 



3/2 



Vk (r) Vfe, 



(139) 



one has the Boltzmann factor 25 j 27V f H>] 7 e q-(117), and 
other physical quantities via new variables ipk- The first 
two terms in eq.(117,) are simple 

/o[V] = 7^73 / (fc z 2 /2 + ar) / ^MwWtt 

= / (fc 2 /2 + a T ) V^Vfc- (140) 
Jk 

The quartic term is 



fint [V>] 



L X Ly 



<J + -0(141) 

7T) JklM'V 



2 (27r)° 

X [k,l|k',l'] 



with 



[k,l|k',l'] 



1 



¥4(r) V f(r)^(r)^(r) (142) 



calculated in Appendix A. Generally the expression is 
not very simple due to the so called " Umklapp" processes 
since when four quasi - momenta involved We turn now to 
the first application of this basis: calculation of harmonic 
excitations spectrum of the vortex lattice. 



For negative ot and neglecting thermal fluctuations 
the minimum of energy is achieved by choosing one of 
the degenerate lattice solutions, the hexagonal lattice <pa 
in our case. This was the main subject of the previous 
section. When thermal fluctuations are weak, one can ex- 
pand in temperature around the mean field solution. The 
zero quasimomentum field is then shifted by the mean 
field solutions. In our new LLL units we therefore ex- 
press the complex fields ^fc via two "shifted" real fields 
(O k = 0*_ k , A k = A*_ k ): 

= v (2tt) 3/2 4 + (O k + iAk) (143) 

with value of the field found in section II in the LLL units 
being 



v o = \—s- 



"Pi" 



(144) 



Notations "O" and "A" indicate an analogy to optical 
and acoustic phonons in atomic crystals. The constants 
Ck will be chosen later and will help to diagonalize the 
quadratic part of the free energy. Substituting this into 
the energy eqs.(140) and (141), one obtains a constant 
"mean field" energy density of section II, 



fmf _ O'T 

~vol ~ ~2/V 
while the quadratic part is 



(145) 



h 



fc 2 /2-a T + 2 V 2 |c k | 2 /3 k (0* k - iA%) 



(O k + iA k ) + ^ / [ 7k c k (0* k + iAi) (Ok + iA k ) + 
4 Jk 



(146) 



c.c 



where functions, 



-l^|^(r)| 2 |^ k (r)| 2 = [0,k|0,k], (147) 

7k = ^^* 2 (r)^ k (r)(^_ k (r) = [0,0|k,-k], 

are calculated and given explicitly in Appendix A. There 
is no linear term, since we shifted by the mean field solu- 
tion. 

The choice 



B. Excitations of the vortex lattice and perturbations 
around it. 

1. Shift of the field and the excitation spectrum 

Shift of the field and diagonalization of the 
quadratic part 



Ck 



1L 

l7k| 



eliminates the OA terms, diagonalizing / 2 : 



1 



h= g j^OlOk+etAlAk. 



(148) 



(149) 
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The resulting spectrum is:epsilon_III 

e° k = /4 k + k 2 j2- = a 2 Ak + fc 2 /2; (150) 
Mok = or + «o ( 2 Pk + l7k|) = (2/3 k + | 7k | - /3 A ) ; 

PA 

/4k = cl t + v 2 (2/3 k - |7k|) = - ^ (2/3 k - | 7k | - Pa) 

Pa 

The cubic and quartic terms describing the anharmonic- 
ities or interactions of the excitations (phonons) are 

h = f S(k z -l z -m z )A 3 (k,l, m )x (151) 

J k,l,m 

[(0* k - iA* k ) (O, + iA t ) (O m + iA m ) + ex.] ■ 
U= I 5(^-/ z + ^-^)A 4 (k,l,k',l') (152) 

Jk,l,k',l' 

x (0£ - iA%) (O, + (Cft - U£,) (Or + *40 , 
where 

A 3 (k, 1, m) = vo 25^/2 t k ' m ] c k c i c m (153) 

A 4 (k, 1, k', 1') = ^ [k, It'll, 1'] (154) 

with [k,k'|l,l'] defined in cq.(142). 

Super soft Goldstone (shear) modes 

While the O mode is "massive" even for small quasi 
- momenta, the A mode is a Goldstone boson resulting 
from spontaneous breaking of several continuous symme- 
tries and is therefore " massless" . The broken symmetries 
include the electric charge U (1) , (magnetic) translations 
and rotations. Spectrum of Goldstone modes is typically 
"soft" and quadratic in momentum. This is indeed the 
case, as far as the field direction z is concerned, eq.(150), 
but the situation in the perpendicular directions is dif- 
ferent (Eilenberger, 1967; Lee and Shenoy, 1972). 

We use expansion of the functions pk and 7 k , eq.(147), 
derived in Appendix A: 

/3k = PA-^k 2 + /?4Ak 4 ; (155) 
7k = Pa - ^yk 2 - ip A k x k y 

+t PAk x k y k* + (k 4 _ 4fc 2 fc 2) 

with constants given in Appendix A, /3 4 a = 0.132. The 

acoustic spectrum consequently has the following expan- 
sion at small momenta: 

l*Ak = ( 2 &A - /?a/8) vl |k| 4 + ... = 0Aa H |k| 4 (156) 

All the quadratic term cancel and the Goldstone bosons 
are " supersoft" . 



One can further investigate the structure of these su- 
persoft modes and identify them with "shear modes" 
(Moore, 1989, 1992; Zhuravlev and Maniv, 1999, 2002). 
To conclude, there are many broken continuous symme- 
tries (translations in two directions, rotations and the 
U (1) phase transformations, forming a rather uncom- 
mon in physics Lie group) leading to a single Goldstone 
mode. The commutators of the magnetic translations 
generators P and the U (1) generator Q = 1 arc (using 
the explicit form eq.(67)): 

[P X ,P V ] = iQ- [P x , Q] = 0; [P y ,Q] - 0, (157) 

and form the so called Hciscnbcrg - Weyl algebra. How- 
ever the Goldstone mode is much softer than the regu- 

i 1 4 i 1 2 

lar one: |k| instead of |k| . The situation is not en- 
tirely unique, since ferromagnetic spin waves, Tkachcnko 
modes in superfluid and excitations in 2D electron gas 
within LLL share this property. A rigorous general 
derivation of the modification of the Goldstone theorem 
in this case is still not available. Note also that, when the 
magnetic part is not neglected, the modes become mas- 
sive via a kind of Anderson - Higgs mechanism, which 
gives them a small "mass" of order 1/k 2 in our units. 

This exceptional "softness" apparently should lead to 
an instability of the vortex lattice against thermal fluctu- 
ations. Indeed naive calculation of the correlator in per- 
turbation theory shows that certain quantities including 
superfluid density \ip\ 2 are infrared (IR) divergent (Maki 
and Takayama, 1971). This was even considered an indi- 
cation that the vortex lattice does not exist (Moore, 1997; 
Nikulov et ai, 1995a; Nikulov, 1995b), despite large body 
of experimental evidence, even at that time. As a result, 
the perturbation theory around the Abrikosov solution 
was not developed beyond the one - loop order for a long 
time. One could argue (Brandt, 1995) that real physics 
is dominated by the small mass 1/k 2 of the shear mode, 
acting as a cutoff that prevents IR divergencies, but basic 
physical properties related to thermal fluctuations near 
H C 2 (T) seemed to be independent of the cutoff, especially 
for high T c superconductors. In (Rosenstein, 1999) the 
IR divergencies were reconsidered and it was found that 
they all cancel exactly at each order in physical quan- 
tities like free energy, magnetization etc. We therefore 
systematically consider the (renormalized) perturbation 
theory for free energy up to two loops and then turn to 
other physical quantities. 

2. Feynman diagrams. Perturbation theory to one loop. 

Feynman diagrams for the loop expansion 

To develop a perturbation theory, the coefficient in 
front of the Boltzmann factor, eq.(117) is considered large 

/ = - [fo + h + h + h] ■ (158) 
a 

The "small parameter" a is actually 1, but will be useful 
to organize the perturbation theory before the actual ex- 
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FIG. 6 Feynman rules for vortex lattice 

pansion parameter is uncovered in the process of assem- 
bling the series. One does not have to consider a linear in 
fields term f\ since it involves only the k = Goldstone 
excitations and does not contribute to bulk energy den- 
sity (Jevicki, 1977). The free energy is calculated from 
eq.(119) by expanding exponent of "vertices" fa [ip] and 
fi [ip], so that all the integrals become gaussian: 

/(aT) = -/ m /- 2 5 / 2 7rlog{ / e -*(A+A+/0} (159) 

= i/ - 2 5 / 2 7rlog{ / e"^[l- I (/ 3 + / 4 ) + ...]} 

= — fo — 2 5 / 2 7rlog{ / e~~-^ 2 } + connected diagrams. 

The propagators entering Feynman diagrams (Fig. 
6a,b) are read from the quadratic part, eq.(140): 



G° > A {k) 



2 5 / 2 7T 
0,A ■ 



(160) 



The leading order propagators are denoted by dashed 
and solid lines for the A and the O modes respectively. 
Nonquadratic parts of the free energy are the three - leg 
and the four - leg "vertices", Fig. 6c-f and Fig. 6g-k 
respectively. It is important for disappearance of "spu- 
rious" IR divergencies (to be discussed later) to realize 
that vertices involving the A field are "soft", namely at 
small momentum they behave like powers of k. For ex- 
ample, the AkAiA m vertex, Fig. 6f, is very "soft". At 
small momenta it is proportional to the fourth power of 
momenta 

The power of a L ~ x , L = \ (3iV 3 + 47V 4 ) - iV 3 - 7V 4 + 
1, where N^,N4 are numbers of the three - leg and the 



four - leg vertices, in front of a contribution means that 
topologically the number of "loops" is L (Itzykson and 
Drouffe, 1991). The leading term, the mean field energy 
is of order a . 

Energy to the one loop order 

Important point to note is that in the "ordered" phase, 
despite the fact that we are talking about perturbation 
theory, the shift or, in other words, definition of the 
"physical" excitation fields Ok and A& in terms of the 
original fields ipk can change from order to order (Itzyk- 
son and Drouffe, 1991). The shift v in eq.(143) is there- 
fore renormalized, that is, 



w 2 + av\ + . 



(161) 



One finds v\ in the same way vq was found, namely, by 
minimizing the effective the free energy at the minimal 
order in which it appears. Let us therefore explicitly 
write several leading contributions to the energy 

/ -/, ' /, ' <>/, ' •••• (162) 

a 

We us start from the "mean field" part in eq.(159): 

o 1 



fmf _ }_ 

vol a 

_ 1 
a 



-axv 



(163) 



1 



a TV + ^0AV O 



+ [-a T v\ + PavIvI] + 



-jAvt 



+ O (a 2 ) . 



The leading order is or 1 and comes solely from the 
mean field contribution, which is therefore the leading 
contribution in eq.(159) and coincides with eq.(145): 



A 

vol 



LLrp 

"vZ' 



(164) 



This part of energy can also be viewed as an equation 
determining vq. 

Substituting vq into the expression in the second 
square bracket in eq.(163) makes it zero. The only con- 
tribution to the order a comes from the second term in 
eq.(159), the "trace log" ,-2 5 / 2 tt log{ / e~^ h } equal to: 



1 



2 3 / 2 7T 2 



{log[G °(fc)] +log[G^(fc)]} (165) 



When we take the leading order in the expansion of the 
excitation spectrum in powers of a 



o,A 



= a T + i; 2 (2/3 k ±|7 k |) (166) 
= a T + v 2 (2/3k ± |7k|) + avj (2/3 k ± | 7k |) + ... 

2 

.2 



(/4k 4 ) +a« 2 (2/3 k ±| 7k | 
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the one loop energy becomes: 

A. = 1 

vol 2 3 / 2 7T 2 



{log[(MSc) 2 + fc 2 2 /2]+log[ 



(4u) 2 + kl/2}} = - [ ( Mo l + Mb1c) = 2.848 \a T \\<&x7) 

If Jk. 



3. Renormalization of the field shift and spurious infrared 
divergencies. 

Energy to two loops. Infrared divergent renor- 
malization of the shift 

To order a, corresponding to two loops, one has the 
first contribution from the mean field part, which con- 
tains V\, namely, the third square bracket in eq.(163). 
The "trace log" term, eq.(165), contributes due to lead- 
ing correction to the excitation spectrum eq.(166): 



1 



23/ V"" 1 
2 1 



2/3 k 



I7k| 



k (Muk) 
2/3 k + |7k| 
M °k 



fcf/2 
Wk 



+ 



2/3k - |7k| 



l7kL 



Mok 



+ *2/2 
(168) 



while the rest of the contributions in eq.(159) are drawn 
as Feynman two - loop diagrams in fig. Fig. 7 and cannot 
contain v\ , since propagators and vertices already provide 
one factor a. The minimization with respect to v 2 results 
in: 



1 



1 



■2/3 k +|7k| , 2/? k -fok| 



2lT Jk ^Ok 



Mo°k 



^Ok 



(169) 



27T/3. 



2/3 k 



I7k| 



M?k 



Due to additional softness of the A mode e^ k oc |k| , 
the first ("bubble") integral diverges logarithmically near 
0: 



d 2 ki 



1 



MokJ 



11/2 



a T (2/?4A//3A-l/8) J 




(f) (g) 



d 2 k , 

2 ^ 1°S ^ 



(170) 

This means apparently that for the infinite infrared cut- 
off fluctuations destroy the inhomogeneous ground state, 
namely the state with lowest energy is a homogeneous 
liquid. It is plausible that since the divergence is loga- 
rithmic, we might be at lower critical dimensionality in 
which an analog of Mermin - Wagner theorem (Itzyk- 
son and Drouffe, 1991; Mermin and Wagner, 1966) is 
applicable. Even this does not necessarily means that 
perturbation theory starting from ordered ground state 
is useless(Jevicki, 1977). A rigorous way to proceed in 
these situations have been found while considering sim- 
pler models like the > 4 model", F = ±(V^ a ) 2 + V(p a2 ), 
in D = 2 with number of components larger then one, say 
a = 1,2. Considering the corresponding statistical sum, 
one first integrates exactly zero modes, existing due to 
spontaneous breaking of a continuous symmetry (U (1) 



FIG. 7 Two loops connected diagrams contributing to free 
energy. 



in our case, field rotations in the ip 4 model) and then de- 
velops a perturbation theory via steepest descent method 
for the rest of the variables. When the zero mode (the 
above mentioned Goldstone boson with k = 0) is taken 
out, there appears a single configuration with lowest en- 
ergy and the steepest descent is well defined. For in- 
variant quantities like energy this procedure simplifies: 
one actually can forget for a moment about integration 
over zero mode and proceed with the calculation, as if 
it is done in the ordered phase. The invariance of the 
quantities ensures that the zero mode integration triv- 
ially factorizes. This is no longer true for noninvariant 
quantities for which the machinery of " collective coordi- 
nates method" should be used (Rajaraman, 1982). In our 
case, we first note that the shift of the field is not a U (1) 
'or translation invariant quantity, so invariant quantities 
like energy might be still calculable. Moreover the sign 
of the divergence is negative and a physically reasonable 
possibility that the shift decreases as a power of cutoff: 



1 — ac log L 



- {ac log L) 2 



(171) 



v n e 



-ac log L 



IR divergences in energy. The "nondiagram- 
matic" mean field and Trlog contributions. 

Substituting the IR divergent correction v\, eq.(169), 
back into the free energy, eqs.(163) and (168), one ob- 
tains a divergent contribution for the " nondiagrammatic" 
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terms in eq. (free_energy_expansion_III) : 
1 /\2/? k 



■pAvj + vf — 

2lT A Mok 

I [, 

2 (2tt) 2 (3 a Vk Mok A 



W + 2/3 k - | 7 k| 



Mok 



(172) 



2/3k 



M? k 



l7k ' " ^^ k )} 2 



other power-wise divergencies left to the two loop order. 
Analogous analysis of the OOA vertex shows that the 
OOA setting sun diagram, Fig. 7c is also convergent. 

Naively logarithmically divergent AAO setting sun di- 
agram, Fig. 7b actually has both the log 2 L and the log L 
divergences. The AAO vertex function is 



containing both the (log L) 

fdiv _ Pa 
vol 



2 (27r) 2 Jk,i M^kMoi 



(173) 



and the sub leading logL divergences. However we 
haven't finished yet with the order a. They also likely to 
have divergences, naively even worse than logarithmic. 
We therefore return to the rest of contributions to the 
two loop order. 

"Setting sun" diagrams. 

One gets several classes of diagrams on Fig. 7, some of 
them IR divergent. The naively most divergent diagram 
fig. 2a actually converges. It contains however two AAA 
vertices, each one of them is proportional to the fourth 
power of momenta. The integrals over k z and l z can be 
explicitly performed using a formula 



1 



1 



1 



1 

27 J kx J h kl/2 + fil PJ2 + (k z + i z y/2 + 

= r- (174) 

MkMlMk+1 (Mk + Ml + Mk+l) 

The divergences appear, when one or more factors in de- 
nominator belong to the A mode for which ,u k oc |k| 2 
for small k. However, if the numerator vanishes at these 
momenta, the diagram is finite. The numerators con- 
tains vertices involving the same " supersoft" field A and 
typically vertices in theories with spontaneous symmetry 
breaking are also soft (this fact is known in field the- 
oretical literature as "soft pions" theorem due to their 
appearance in particle physics). In the present case they 
are "super soft". The AAA vertex function, eq.(154), is 



, AAA 



(k,l,m) 



iL x L y v 



Rc lCm [0,-k|l,m] (175) 



A AAO (k, 1, m) = V -^{ ClC * m ci [-k, -HO, 1] (177) 

+4 c i c m [0, -k|l, m] + c,c k c m [0, -l|k, m] - <„cic k x 
[0, -m|l, k] - c,c m c k [-1, -k|0, m] + cic^c k [-1, -m|0, k]}. 

We will need its asymptotic when one of the momenta of 
the soft excitation A is small 



A AAO (0,1, m) 



2 2 n 3/2 



l7i|- 



(178) 



The diagram of Fig. 7b, after integration over the field 
direction momenta k Zl l Zl is: 



2VL 2 



A AAO (k, 1, m) A AAO (-k, -1, -m) 
kj,m Mk MlVm (M k + li-x + M? ^ 



(179) 

The leading divergence is determined by the asymp- 
totics of the integrand as both k and 1 approach zero. 
Consequently it is given by the integral when the two 
vertex functions replaced with their values taken at k = 
1 = and momenta of u£ and (/z k + u A + a^) in the 
denominator also taken to zero. The log 2 L divergent 
part near k = 1 = is therefore 



fdiv = -2 5 7r 3 i 2 



A AAO (0, 1, m) A AAO (0, -1, -m) 
k,i,m MkMiVP (Mo + Mi 4 + Mp) 



L> z L x Ly 2 



2|7o| 



Q2_2 / A..A..O..O 

1 Jk,i M k Mi Mo Mo 

vol f Pa 
(2tt) 2 7 k i 



M k Mf 



(180) 



2 5 7ri 3 

+c 1 *c k c m [0, -l|k, m] + c^cic k [0, -m|l, k] - c k c, c^x 
[-1, -m|0, k] - qc k c^ [-k, -m|0, 1] - c m c, c k [-1, -k|0, m]} 

One easily sees that for each of the " dangerous" momenta 
k = 0, 1 = orm = each one of two vertices vanishes. 
For example when k = 



The "bubble" diagrams and cancellation of the 
leading divergences 

Diagrams given in Figs .7e,f,g, can be easily evaluated: 



/(., 



f,g) 



vol (2tt) 2 J k ,i 



L x L y v 



A AAA (0, 1, m) = * 25^7/2 3 



^{cic m [0,0|l,m] 



c*c m [0, -1|0, m] + c* mCl [0, -m|l, 0] - c\c* m [-1, -m|0, 0] 
-ciC [0, -m|0, 1] - c m c\ [-1, O|0, m]} (176) 

This means that there are at least two powers of k in 
the numerator and the integral converges. There are no 



Wa (2tt 
The leading divergence is 
3 1 



— -M/ |7kl u-i 



Jdiv 
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vol 2 (2tt) 2 ./k/'^MkMi 4 ' 



(181) 



(182) 



One observes that sum of three leading (logL) 2 diver- 
gences given in eqs. (173), (180) and (182) cancel. There 
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are still sub leading log L divergences. They require more 

care, since "not dangerous" momenta cannot be put to 

zero, and are treated next. 

Cancellation of the IR divergencies 

The two - loop contribution to energy in a "standard" 

form: 



/= 



V 



(27T) 2 i k ,l msX 



(183) 



In order to demonstrate cancelation of the IR diver- 
gences we investigate the value of the numerator F (k, 1) 
at k = and 1 = and show that F (k = 0, 1) = 
and F (k, 1 = 0) = 0. The part due to nondiagrammatic 
terms eq.(172) can be written as: 



F 1 (k,l) 



1 r 2/3 k +|7k| 



Mp 



Mk + 2/?k ~ |7k|] 



r 2/3i+l7i 
[ Mp 



Hi + 2/?i 



F > (0 , 1) = _^[?^N 

similarly, the setting sun diagram, 
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(184) 
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(185) 



according to eq.(180). The divergent part of the bubble 
diagrams can be written as: 

F 3 (0, 1) = ^[| 7 ,| + 2/3,(^4 - (186) 



Mf Mf 



Mf Mf 



One explicitly observes that F 1 (0, 1) + F 2 (0, 1) + 
F 3 (0, 1) = . The same happens F 1 (k, 0) + F 2 (k, 0) + 
F 3 (k, 0) = . Therefore all the IR divergences, e.g. the 
logi and (logL) , cancelled. Similar cancellations of 
all the logarithmic IR divergencies occur in scalar mod- 
els with Goldstone bosons in D = 2 and D — 3 (where 
the divergencies are known as " spurious" (David, 1981; 
Jevicki, 1977)). 

Vortex lattice energy 

The finite result for the Gibbs free energy to two loops 
(finite parts of the integrals were calculated numerically) . 
Up to two loops the calculation(Li and Rosenstein, 
2002a, b,c) (extending the one carried in rcf. (Rosenstein, 
1999) to Umklapp processes) gives: 



fd (or) = 
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2.848 \a T 



1/2 



vol 2f3 A 
In regular units the free energy density is 
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(187) 
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2.848 \a T \ 1/2 + — 
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(188) 



Below we use this expression to determine the melting 
line and various thermal and magnetic properties of the 
vortex solid: magnetization, entropy, specific heat. Near 
the melting point ar ~ —9.5 the precision becomes of 
order 0.1% allowing comparison with the free energy of 
vortex liquid, which is much harder to get. Eventually 
the (asymptotic) expansion becomes inapplicable near 
the spinodal point at which the crystal is unstable due 
to thermal " softening" . This is considered using gaussian 
approximation in subsection D. 



4. Correlators of the U (1) phase and the structure function 

Correlator of the U (1) phase and helicity mod- 
ulus 

The correlator of the order parameter is finite at all 
distances in the absence of thermal fluctuations 

C mf (r, /) = r (r) V (/) = « V (r) <p (/) , (189) 

where v 2 , = jj^- is finite exhibiting the phase long range 
order in the vortex lattice (despite periodic modulation). 
However the order is expected to be disturbed by thermal 
fluctuations. Leading order perturbation theory gives an 
early indication of the loss of the order in directions per- 
pendicular to the field. The leading correction consists of 
the a correction to the shift v, eq.(169) and sum of two 
propagators 

C(r,r')= f V * k {r)^{r')x (190) 

Jk,l 

<l * + ^^ ( °'- ! ' 4l,1W + ^P 

x (Oj + iAi)] >= (r, r') + a{v 2 lV * (r) <p (/) + 

— J-3 / [Gq (k) + G„ (k)] <p% (r) Vk (/)} + O (a 2 ) . 
Z{2ir) Jk 

One observes that the logarithmic divergences of the sec- 
ond and the third terms cancel, but that the correlator 
in the x — y plane depends on the large distance r — r' as 
a log: 



C(r,z = 0,r',z' = 0)~i; o V W^(r') 

; ik-(r-r') _ ^ ' 



(191) 



Mifk } 



= vltp* (r) tp (r') [1 + aclog |r - r'|] . 

It is expected, that exactly as the expectation value v 
dependence on IR cutoff, the actual correlator is not 
growing logarithmically, by rather decreasing as a power 
|r — r'p c . This is an example of the Berezinsky - Koster- 
litz - Thouless phenomenon (Itzykson and Drouffe, 1991). 
It appears however at rather high dimensionality D = 3. 
Note however that the LLL constraint (large magnetic 
fields) effectively reduce dimensionality, enhancing the 
role of thermal fluctuations. 
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In the direction parallel to the field the correlations are 
still long range. Indeed the helicity modulus is 



C(r = 0,z,r' = 0,z') 



(192) 



{1 + 



2(2nf 



D ikz (z — z') 



} 



Structure function. Definitions 

The superfluid density correlator, defined by 



S{r) 



= -[ 

vol J r , 



\4>{r')\ 2 \i>{r' + r)\' 



(193) 



quantifies spontaneous breaking of the translational and 
rotational symmetries only as in both locations the su- 
perfluid density is invariant under the U (1) gauge trans- 
formations. This is different from the phase correlator 
(ip* (r) ip (/)) discussed in the previous subsection, which 
decays as a power as indicated by the IR divergences. As 
in the case of the U (1) phase correlations, it is easier 
to consider the Fourier transform of the correlator, the 
structure function. Since translational symmetry is not 
broken along the field direction, one can restrict the dis- 
cussion to the lateral correlations and consider partial 
Fourier transform: 



5(q,0) - j dre^ r 5(r, r z = 0) 



(194) 



In this subsection the structure function is calculated 
to leading order in thermal fluctuations strength (har- 
monic approximation) within the LLL, namely neglect- 
ing higher an corrections. We discuss these corrections 
later. 

Structure function of the vortex crystal without 
thermal fluctuations. 

Substituting the LLL mean field solution cq.(61) into 
the definition of structure function one obtains: 
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K 

where we made use of formulas in Appendix A and delta 
function peaks are located at vectors of the reciprocal 
lattice. The height of the peak decreases rapidly be- 
yond the reciprocal magnetic length (which is our unit). 
When mesoscopic thermal fluctuations are significant, 
they might broaden the peaks far below the tempera- 
ture at which the lattice becomes unstable (the spinodal 
point). 

Leading order corrections to thermal broaden- 
ing of Bragg peaks 

The calculation of the structure function closely follows 
that of free energy. The correlator is calculated using the 



Wick contractions: 



5(r,,z = 0) 



L X Ly 



y J k,l,k' ,1' ,r' 



xip r (r' + r)[v6 k + 



V2(2tt) 



^(r>,(r>£,(r' + r) 

(O k - iA k )}[vSi 



3/2 



(O k > - iAy)][v5 v + 



Cl' 



V2(2tt) 



3/2 



(O v +iA v )]. (196) 



The leading order (a ) term is the mean field part, 
eq.(195), while the first order term is the harmonic fluc- 
tuation part. 

The fluctuation part contains v\ corrections term 
54 (q, 0), same in structure as the leading order but 

with (IR diverging) coefficient 2 vfAn 2 instead of 
(jf^j 47r 2 and four contractions (diagrams): 



^(r>,(r>£,(r'+r) 



L x L y [V2(2irf /2 } 2 Jk,i,k',i',r' 
x^(r' + r)K 2 (G& - G$ k ) SfoSw + (197) 
c 2 i (Got - Got) SkSk'Si+v + 2 (G ° fe + G&) 5 k '5 v 5 k+l 
+2{G$ k + G$ k )5 l & k .8 k+v }. 

Performing integrations and Fourier transforms using 
methods described in Appendix A the first two contri- 
bution are: 



5i(q,0) = 



xe 2 



cos (fc x fc„ + k x Q + 6> k ) (198) 
Wiir 1 - GO"' 



where Q and k are the "integer" and the "fractional" 
parts of q, in a sense q = k+ Q for k inside Brillouin 
zone and Q on the reciprocal lattice. The third term is 



5 2 (q,0) 



AirciT 



while the last is: 



5 3 (q,0) 



^n(q)e" 



cos(k x Q) (200) 



W + i^y 1 



Cancellation of the infrared divergences 

Although all of the four terms 5i, 52, S3 and 54 are 
divergent as any of the peaks is approached, k — > 0, the 
sums 5i,52 and 53,54 are not. We start with the first 
two: 



5i+5 2 = -5— e 2 /i(q) 

PA 



(201) 



26 



where 

/i(q) = [1 + cos (k x k y + k x Q + c fe )] (^y 1 
+ [1 - cos(k x k v + k x Q + Ck)] (^ok) _1 (202) 

When k — > 0, it can be shown that k x k y + Ck = O (k 2 ' 2 ) 
, thus (k x k y + kxQ + Cfe^kxQ, and 1 — cos(k x k y + 
k x Q + Cfe) — > (k x Q) 2 . Hence it will cancel the 1/fc 2 
singularity coming from l/VoV Thus /i(q) approaches 
const. + const. ■ ^ k> ^ when Q ^ 0, and approaches 
const. + const. ■ k e when Q = 0. Similarly the sum 
of 54(q,0) and 53 (q, 0) is not divergent, although sep- 
arately they are. Their sum is. 



4-kci 



1/2 



-^„(q) e -^[/ 2 (Q) + / 3 ], 

PA 

(203) 



5 3 (q,0) + 5 4 (q,0) = 
with 

/ 2 (Q) = ^[-l + cos(kxQ)][( Mo ° k )" 1 + (M4) _1 



h 



/ (MoL + Mok) 



-8.96 



(204) 



Super so ft phonons and the "halo" shape of the 
Bragg peaks 

The sum of all the four terms can be cast in the fol- 
lowing form: 



S(q,0) = -^-f<S n (q)e-V 

PA PA 



(205) 



[/i(q) + 5„(q)/ 2 (Q) + 5 n (q)/s] 5 



The results were compared(Li and Rosenstein, 2002a) 
with numerical simulation of the LLL system in (Sasik 
and Stroud, 1995). For reciprocal lattice vectors close to 
origin the value of /2(Q) are: 

Table 1 

Values of / 2 (Q) with small ni ; ri2. 



ni,n 2 


(0,1), (1,0), (1,1) 


(1,2), (0,2), (2, 2) 


1,3 


/ 2 (Q)/(2tt) 


-5.20 


-7.11 


-8.31 



The correction to the height of the peak at Q, 
(l+^M / 2 (Q)> i s quite small. The theoretical prediction 



(l+ci/s)- 

has roughly the same characteristic saddle shape " halos" 
around the peaks as in MC simulation ref. (Sasik and 
Stroud, 1995) and experiment (Kim et at, 1999). Con- 
versely, MC simulation result provides the nonperturba- 
tive evidence u A k — + |fc| for small k x ,k y . In eq.(205), 
if Mok — * 1^1 7 we would get a contribution from the 



most singular term const. + const. 



(kxQ) 



This term 



will become constant when k — ► 0, and we will not get 
the same characteristic saddle shape "halos" around the 
peaks as in ref. (Sasik and Stroud, 1995). Consequently 
the /ig k — > |fc| asymptotics for k — > is crucial for such 



characteristic shape and thus the MC simulation result 
provides a nonpcrturbativc evidence for it. 
Magnetization profile 

Another quantity which can be measured is the mag- 
netic field distribution. In addition to constant magnetic 
field background there are 1/k 2 magnetization correc- 
tions due to field produced by supercurrent. To leading 
order in 1/k 2 magnetization is given by eq.(123). The 

2\ • 



superfluid density (^\'ip(r)\ j is calculated as in eq.(190): 

|V(r, z = 0)| 2 ) = f- Mr)\ 2 + -L / |^kW| 2 (4- + 4r) 
/ (3 A 2tt 7 k ^ ^ k 



J_ /' f 2/3k + |7kl 2/? fc -j 7fc U , , , 



'> Mok Mok 

i|^(.r)|v (206) 



Its Fourier transform p(q) = J dre iqr ^|^>(r, z = 0)|' 
can be easily calculated: 



p(q) = 4. 2 <5 n (q){^ + T l /[ (e *x q _2A + W 



'Ai 2^ 



ft 



) 



x(^)"V(e^-^M)( Alo 4 k) - 1 ]} 

Pa 



xe 



-^+i-7mi(n 2 + l) 



(207) 



Performing integrals, one obtains: 

p(q) = 47r 2 (5„(q)e-T+^™i(™2+i) x p og j 

{g + ^[/3 + MQ)]^ 1/2 } 

The function / 2 (Q) and constant /3 appeared in eq.(204). 



C. Basic properties of the vortex liquid. Gaussian 
approximation. 

1. The high temperature perturbation theory and its 
shortcomings 

The loop expansion 

Unlike the perturbation theory in the crystalline state, 
in which various translational, rotational and gauge U (1) 
symmetries are spontaneously broken, the perturbation 
theory at high temperature is quite straightforward. One 
directly uses the quadratic and the quartic terms in the 
Boltzmann factor, eq.(117) as a "large" part K and a 
"perturbation" V : 



K = 



2 5 / 2 TT 



2 5 / 2 7T' 



(209) 



Again the "parameter" a is actually 1, but is regarded 
as small and the actual expansion parameter will be- 
come apparent shortly. The Feynman rules for a field 
ip, namely the propagator 



G (k) 



2 5 / 2 7T 

fc 2 /2 + a T 



(210) 
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FIG. 8 (a),(b)Feynman rules in the homogeneous phase;(c), 
the two loop correction to energy 

and the four - point vertex are given on Fig. 8a, b re- 
spectively is defined in eq.(142). Since ip is a complex 
field, we use an arrow to indicate the "orientation" of 
the propagator. 

The leading contribution to the LLL scaled free energy, 
that of the quadratic theory, is 

A = 2 5 /2 7rT rlog J D- 1 (fc) = ^^^logG(fc) (211) 
vol f , fc 2 /2 + a T , . i/2 

= ^172-2 l0 S 2 5/ 27T = V ° l 4a T + ™™t 

The const in the last line is ultraviolet divergent, but 
unimportant for our purposes and will be generally sup- 
pressed. Corrections can be conveniently presented as 
Feynman diagrams. 

The next, "two - loop", correction is the diagram in 
fig. 8c, which reads 

^ = tAf / A(M,M)G o (fc)Go(0 (212) 
(2n) Jk,l 

= vol Aaj} 

Observe that the k z integrations can be reduced to cor- 
responding integrations in quantum mechanics of the an- 
harmonic oscillator (Ruggeri and Thouless, 1976; Rug- 
gcri, 1978; Thouless, 1975), so that the series resemble a 
dimensionally reduced D — 2 = l field theory or quantum 
mechanics. 

Actual expansion parameter and the applicabil- 
ity range 

This expansion can be carried to a high order after sev- 
eral simple tricks are learned (Hu and MacDonald, 1993; 
Hu et al, 1994; Ruggeri and Thouless, 1976; Ruggeri, 
1978; Thouless, 1975). The result to four loops is: 

^ = 44 /2 + -- J i + S- (213) 
«t 2a?/ 2 24 4 

One observes that the small parameter is a T 3 ^ 2 although 
coefficients grow and series are asymptotic. The differ- 
ence with analogous expansion in the crystalline phase 



is that the sign of cit is opposite and the leading or- 
der is t/wt rather than a\. Phenomenologically the re- 
gion of positive large ar is not very interesting since at 
that point, for example, magnetization is already very 
small. Also higher Landau levels effects become signif- 
icant as will be discussed in subsection E, where HLL 
effects (Prange, 1969) are considered. 

Therefore attempts were made to extend the series to 
smaller temperatures. One of the simplest methods is to 
perform a Hartree - Fock type resummation order by or- 
der. Let us first describe in some detail a certain variant 
of this type of approximation called generally gaussian, 
since it will be extensively used to treat thermal fluctua- 
tions as well as disorder effects in the following sections. 
It will be shown in subsection D that the approximation, 
is not just a variational scheme, but constitutes a first ap- 
proximant in a convergent series of approximants (which 
however are not series in an external parameter) called 
"optimized perturbation theory" (OPT). 

2. General gaussian approximation 
Variational principle 

We start from the simplest, one parameter version of 
the gaussian approximation which is quite sufficient to 
describe the basic properties of the vortex liquid well 
below the mean field transition point ax — 0. Within 
this approximation one introduces a variational param- 
eter u (which is physically an excitation energy of the 
vortex liquid) adding and subtracting a simple quadratic 
expression /i 2 |i/>| 2 from the Boltzmann factor: 

f{fj)=K + aV (214) 
K = ^/(^H 2 + ^| 2 ) = jUiEG-Vfc(215) 

v = _±-[a T ^i^i 2 + i J k A(k,i,k',i')r k TpirM 

where the constant a was defined by 

a = a T -n 2 . (216) 

Now one considers K as an "unperturbed" part and aV 
as a small perturbation. This is a different partition than 
the one we used previously to develop a perturbation the- 
ory Despite the fact that a = 1, we develop perturbation 
theory as before. To first order in a, the scaled free en- 
ergy is: 

W = -2 5/2 7 rlog{^ e -^} (217) 

= -2 5 ^ log{ f e- K+aV } w -2 5 ^ log{ [ e- K [1 + aV}} 

Jifi Jtf> k 

= -2 5 ^tt log{Z M + a [ e- K V} » -2 5 / 2 7r[log Z^-a (V)^] 

J 4>k 
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A" 



where Z M is the gaussian partition function Z a = ^ 
and thermal averages denoted by (..) are made in this 
quadratic theory. 

Collecting terms, one obtains 



gauss 

vol 



2fi + a 



2/j, 



2a 4 
— + — 



(218) 



Now comes the improvement. One optimizes the solvable 
quadratic large part by minimizing energy for a = 1 with 
respect to fi. The optimization condition is called "gap 
equation" , 



0. 



(219) 



since the BCS approximation is one of the famous appli- 
cations of the general gaussian approximation. 

Existence of a metastable homogeneous state 
down to zero temperature. Pseudocritical fixed 
point. 

It is clear that the overheated solid becomes unstable 
at some finite temperature. It not clear however whether 
over - cooled liquid becomes unstable at some finite tem- 
perature (like water) or exists all the way down to T = 
as a metastable state. It was shown by variety of meth- 
ods that liquid (gas) phase of the classical one compo- 
nent Coulomb plasma exists as a metastable state down 
to zero fluctuation temperature with energy gradually 
approaching that of the Madelung solid and excitation 
energy diminishing (Leote de Carvalho et at, 1999). It 
seems plausible that the same would happen with any 
system of particles repelling each other with sufficiently 
long range forces. In fact the vortex system in the Lon- 
don approximation becomes a sort of repelling particles 
with the force even more long range than Coulombic. 

Note that there always exists one solution of this cu- 
bic equation for positive [i for all values of ax, negative 
as well as positive. The excitation energy in the liquid 
decreases asymptotically as 



Ma 



4 

ax 



(220) 



at temperatures approaching zero. Importantly it be- 
comes small at the melting point located at ax — —9.5, 
see below. The gaussian energy is plotted on Fig. 9 
(marked as the TO line). The existence of the solution 
means that the homogeneous phase exists all the way 
down to T = albeit as a metastable state below the 
melting point at which the free energy of the solid is 
smaller, see Fig. 11. Physically this rather surprising 
fact is intimately related to repulsion of the vortex lines. 
It is well known that if in addition to repulsion there 
exists an attraction like a long range attractive forces 
between atoms and molecules, they will lead to a spin- 
odal point of the liquid (Lovett, 1977). However, if the 
attractive part is absent like in, for instance, electron liq- 
uid, one component plasma etc., the spinodal point is 
pushed down to zero temperature. It becomes a "pseud- 
ocritical" point, namely, exhibits criticality, but globally 
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FIG. 9 Free energy in liquid. The curve TO is the gaussian 
approximation, while Tl,... are higher order renormalized 
perturbation theory results. Optimized perturbation theory 
gives curves 1, 2, ... and finally BP lines are the Borel - Pade 
results. 



unstable due to existence of a lower energy state (Com- 
pagner, 1974). Scaled LLL free energy density diverges 
as a power as well 



-oo)~- T , 



(221) 



see Fig. 9. 

Assuming absence of singularities on the liquid branch 
allows to develop an essentially precise theory of the LLL 
GL model in vortex liquid (even including overcooled liq- 
uid) using the Borel - Pade (BP) (Baker, 1990) method 
at any temperature. This calculation is carried out in 
subsection D4. The gaussian liquid state can be used 
as a starting point of "renormalized" perturbation the- 
ory around it. Such an expansion was first developed 
by Ruggeri and Thouless (Ruggeri and Thouless, 1976; 
Ruggeri, 1978; Thouless, 1975) for the GL model. 



D. More sophisticated theories of vortex liquid. 

1. Perturbation theory around the gaussian state 

After the variational spectrum fi was fixed, one can ex- 
pand in presumably small terms in eq.(214) multiplied by 
a up to a certain order. Here we summarize the Feynman 
rules. 

Feynman diagrams 

The propagator in the quasi - momentum space, 



Gfc = 



M 2 + fc 2 2 /2' 



(222) 



is the same as in usual perturbation theory, Fig. 8a, but 
with gaussian mass fi. The four - leg interaction vertex 
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is also the same as Fig. 8b, but there is an additional 
two - leg term. It has a factor a and treated as a vertex 
and can be represented by a dot on a line, Fig. 10a, with 
a value of 25 " 27r a. The second term is a four line vertex, 
Fig. 8b, with a value of 27 %^ - 

To calculate the effective energy density / = 
— 2 5 / 2 7rlnZ, one draws all the connected vacuum dia- 
grams. To the three loop order one has: 

h = -^ 17+8fl " +aV ) ; (223) 

= — !-r (907 + 510a^ + 96aV + 6a V) ■ 
vol 24^i 8 v 

The liquid LLL (scaled) free energy is generally written 
as (using the gap equation) 

-L=^[l+g(x)\. (224) 
The function <? can be expanded as 

where the high temperature small parameter is x = 
i/i~ 3 . The coefficients c„ which were calculated to 6 t?l 
order in (Ruggeri and Thouless, 1976; Ruggeri, 1978; 
Thouless, 1975) and extended to 9 th order in (Brczin 
et al, 1990; Hikami et al, 1991). The consecutive ap- 
proximants are plotted on fig. 6 (Tl to T9). 

Applicability range and ways to improve it 
One clearly sees that the series are asymptotic and 
can be used only at ar > —2. Therefore the great effort 
invested in these high order evaluations still falls short 
of a required values to describe the melting of the vor- 
tex lattice. One can improve on this by optimizing the 
variational parameter /i at each order instead of fixing it 
at the first order calculation. This will lead in the fol- 
lowing subsections to a convergent series instead of the 
asymptotic one. The radius of convergence happens to 
be around ar — — 5 short of melting and roughly at the 
spinodal point of the vortex solid (see next subsection). 

Another direction is to capitalize on the "pseudocrit- 
ical fixed point" at zero temperature. Indeed, the exci- 
tation energy, for example, behaves as a power u oc a^ 1 , 
other physical quantities are also "critical", at least ac- 
cording to gaussian approximation. It is therefore pos- 
sible to consider supercooled liquid or liquid above the 
melting line but at low enough temperature as being in 
the neighborhood of a pseudocritical point. To this end 
the experience with critical phenomena is helpful. One 
generally develops an expansion around a weak coupling 
unstable fixed point (high temperature in our case) and 
" flows" towards a strong coupling stable fixed point (zero 
temperature in our case) (Itzykson and Drouffc, 1991). 
Practically, when higher order expansions are involved, 
one makes use of the renormalization group methods in 
a form of the Pade - Borel resummation (Baker, 1990). 



This route will be followed in subsection 3 and will lead to 
a theory valid for arbitrarily low temperature. The OPE 
will serve as a consistency check on the upper range of 
applicability of the resummation, which is generally hard 
to predict. 



2. Optimized perturbation theory. 

General idea of the optimized gaussian pertur- 
bation theory 

We will use a variant of OPT, the optimized gaussian 
series (Kleinert, 1995) to study the vortex liquid(Li and 
Rosenstein, 2001, 2002a, b,c). It is based on the "principle 
of minimal sensitivity" idea (Okopinska, 1987; Stevenson, 
1981), first introduced in quantum mechanics. Any per- 
turbation theory starts from dividing the Hamiltonian 
into a solvable " large" part and a perturbation. Since we 
can solve any quadratic Hamiltonian we have a freedom 
to choose "the best" such quadratic part. Quite generally 
such an optimization converts an asymptotic series into 
a convergent one (see a comprehensive discussion, refer- 
ences and a proof in (Kleinert, 1995)). The free energy is 
divided into the "large" quadratic part and a perturba- 
tion introducing variational parameter u like for gaussian 
approximation, eq.(214), although now the minimization 
will be made on orders of a higher than the first. 

Expanding the logarithm of the statistical sum to order 
a n+1 

Z= f cxp(-K) cxp(-aV) = { J2^( aV y e M-K), 
7 n M = -2 5 / 2 7rlogZ = -2 5 / 2 7r (226) 
x{log[/ c -*] - £ L£L (v*)^, 

where () K denotes the sum of all the connected Feynman 
diagrams with G as a propagator and then taking a — ► 1, 
we obtain a functional of G. To define the n th order OPT 
approximant /„ one minimizes f n [G] with respect to G: 

/„ - min/>]. (227) 

M 

Till now the method has been applied and comprehen- 
sively investigated in quantum mechanics only ((Kleinert, 

1995) and references therein) although attempts in field 
theory have been made (Belief et al, 1996a, b; Bender 
et al, 1994; Duncan and Jones, 1993; Guida et al, 1995, 

1996) . 

Implementation and the convergence radius in 
GL 

We can obtain all the OPT diagrams which do not ap- 
pear in the gaussian theory by insertions of bubbles and 
the additional vertex figlc. insertions from the diagrams 
contributing to the non - optimized theory. Bubbles or 
"cacti" diagrams, see Fig. 8, are effectively inserted into 
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FIG. 10 Additional diagram of the renormalized perturba- 
tion theory shown in (a). Bubbles or cacti diagrams summed 
by the optimized expansion are shown in (b) - (d). A diagram 
which is not of that type is shown in (f). 



energy by a technique known in field theory (Kleinert, 
1995). One writes / in the following form: 



f=4fj, 1 +Afj, 1 f(x), 



(228) 



where x = —372 and \i\ is given by a solution of cubic 
equation 



Mi — M2M1 — 4a = 0. 



(229) 



Summing up all the insertions of the mass vertex, which 
now has a value of 25 %^ a, is achieved by 



e 2 = £ + aa. 



(230) 



We then expand / to order a n+1 , and then taking a = 1, 
to obtain f n . The solution of eq.(229) can be obtained 
perturbatively in a: 



Hi = Hi + —T - 



2a 6a 2 32a 3 210a 4 



H2 Hi 



+ 



M2 



M2 1 



+ 



(231) 



The nth OPT approximant f n is obtained by minimiza- 
tion of f n (n) with respect to fi: 



8 



d_ 

da 



fn (M) a ) = °- 



(232) 



The above equation is equal to ^j~( 3 ™+ 4 ) times a polyno- 
mial g n (z) of order n in z = afj,. That eq.(232) is of this 
type can be seen by noting that the function / depends 
on combination . , , - — only. We were unable to prove 

this rigorously, but have checked it to the AQ th order in a. 
This property simplifies greatly the task: one has to find 
roots of polynomials rather than solving transcendental 
equations. There are n (real or complex) solutions for 
g n (z) = 0. However (as in the case of anharmonic oscil- 
lator (Kleinert, 1995)) the best root is the real root with 
the smallest absolute value,. 



We then obtain /U solving the cubic equation, z n — 
a/j, = {clt — /i 2 ) /Lt, explicitly: 

H = 2 1/3 a T ^-27z+ ^-1084 + 729^ (233) 

+ 32^ (~ 2?Z + V - 108a T + 729z2 ) 1 ■ 

For zo — —4, we obtain the gaussian result, dashed line 
marked "TO" on Fig. 9. 

Feynman rules undergo minor modifications. The mass 
insertion vertex, now has a value of ag, while the 



four line vertex is 



25/2 T 



2 S/2 7 



However since the propagator in 



the field direction z and perpendicular factorizes, the k z 
integrations can be reduced to corresponding integrations 
in quantum mechanics of the anharmonic oscillator, as 
we explained in subsection B. Expanding / in a to order 
n+ 1, then one then sets a = 1 to obtain f n . We list here 
first few OPT approximants f n : 



4^ 



2a H 



fo 

fi = (l7 + 8a H [i + a 2 Hl i 2 ); 

h 



M 
1 

2~^ 
1 



(234) 



h + (907 + 510a H ^ + 96a 2 H n 2 + 6a%n 3 ) , 



with higher orders given in ref.(Li and Rosenstein, 
2002b). 

Rate of convergence of OPT 

The remarkable convergence of OPE in simple mod- 
els was investigated in numerous works (Bellet et al., 
1996a, b; Bender et al., 1994; Duncan and Jones, 1993; 
Guida et al., 1995, 1996). It was found that at high 
orders the convergence of partition function of simple in- 
tegrals (similar to the "zero dimensional GL" studied in 
(Wilkin and Moore, 1993) ), 



-f 

J — c 



dipe' 



(235) 



is exponentially fast. The reminder in bound by (Bellet 
et al., 1996a, b; Bender et al., 1994; Duncan and Jones, 
1993; Guida et al, 1995, 1996) 



rjv = \Z — Z N \ < a exp[— C2-2V]. 



(236) 



For quantum mechanical anharmonic oscillator (both 
positive and negative quadratic term) it is just a bit 
slower: 



R 



N 



\E-E N \ <c 1 cxp[- C2 iV 1 / 3 ], 



(237) 



where E is the ground state energy. We follow here the 
convergence proof of (Bellet et al., 1996a, b; Bender et al, 
1994; Duncan and Jones, 1993; Guida et al, 1995, 1996). 
The basic idea is to construct a conformal map from the 
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original coupling g to a coupling of bounded range and 
isolate a nonanalytic prcfactor. Suppose we have a per- 
turbative expansion (usually asymptotic, sometimes non 
Borel summable) 



E(9) = £ c n g n . 



(238) 



n=U 



One defines a set of conformal maps dependent on pa- 
rameter p of coupling g onto new coupling (3 : 



(239) 



While range of g is the cut complex plane the range of /3 is 
compact. The value of parameter p for each approximant 
will be defined later. Then one defines a "scaled" energy 



n(f3,p) = (l-pyE(g([3,p)), 



(240) 



where the prefactor (1 — [3) a is determined by strong 
coupling limit so that Q(f3,p) is bounded everywhere. 
Approximants to Q are expansion to iVth order in /3, 



N 



1 d 



a N (Jt,p) = Y,n\W [{1 ~ ^ E ^M , (241) 



n=0 



with parameter p substituted by p = -|(1 — /3) K . The 
energy approximant becomes 



E N {fl) 



(242) 



Two exponents a = \ an d ft = § , for example, an- 
harmonic oscillator and 3D GL model. OPE is equiv- 
alent to choosing (3 which minimizes E N (f3). It can be 
shown quite generally (see Appendix C of paper in (Ben- 
der et al., 1994) and (Kleincrt, 1995)) that the minimiza- 
tion equation is a polynomial one in p . This is in line 
with our observation in the previous subsection that min- 
imization equations are polynomial in z with p identified 
as —j. 

The remainder Rn = \E — E^\ using dispersion rela- 
tion is bounded by 



Rn < c ig °/*(pNT + c 2 cxp[-N [ i- 



b\N 



1/k 



(243) 



where exponent b is determined by discontinuity of E(g) 
at small negative g: 



Disc E(g) ~ exp[— 



const . 



(244) 



The constants are b = 1 for anharmonic oscillator and 
6 = 3/4 for 3D GL model (Ruggeri and Thouless, 1976; 
Ruggeri, 1978; Thouless, 1975). For 3D GL model, we 
found that Rn < c\ exp (— c 2 A r1 / 3 ) , as in anharmonic 
oscillator. 



3. Overcooled liquid and the Borel - Pade interpolation 
Borel - Pade resummation 

We have already observed using the gaussian approx- 
imation that there exists a pseudo - critical fixed point 
at zero fluctuation temperature c<t — * — oo. One can 
therefore attempt to use the RG "flow" from the weak 
coupling point, the perturbation at high temperature to 
this strongly couple fixed point. This procedure always 
have an element of interpolation. It should be consis- 
tent with the perturbation theory, but goes far beyond 
it. Technically it is achieved by the Borel - Pade (BP) ap- 
proximants. We will not attempt to describe the method 
in detail, see textbooks (Baker, 1990), and concentrate 
on application. 

The procedure is not unique. One starts from the 
renormalized perturbation series of g(x), calculated in 
subsection B, eq.(225), g (x) = ^c n x n . We will denote 
by gk (x) the [k, k — 1] BP transform of g(x) (other BP 
approximants clearly violate the correct low temperature 
asymptotics and are not considered). The BP transform 
is defined as 



f 

Jo 



g' k (x t) exp (— t) dt 



(245) 



where g' k is the [k, k — 1] Pade transform of the better 
convergent series 



(246) 



71=1 



The [k, k — 1] Pade transform of a function is defined as 
a rational function of the form 



(247) 



whose expansion up to order 2fc — 1 coincides with that 
of the function the series eq.(246). 

The results are plotted on Fig. 9 as solid lines for 
k = 3, 4 and 5. The lines for k = 4, 5 are practically indis- 
tinguishable on the plot. The energy converges therefore, 
even at low temperatures below melting. It describes 
therefore the metastable liquid up to zero temperature. 
Due to inherent non - unique choice of the BP approxi- 
mants it is crucial to compare the results with convergent 
series (within the range of convergence). This is achieved 
by comparison with the OPT results of the previous sub- 
section. 

Comparison with other results 

As is shown on Fig. 9, the two highest available BP 
approximants are consistent with the converging OPT 
series described above practically in the whole range of 
a T . One can compare the results with existing (not very 
extensive) Monte Carlo simulation and agreement is well 
within the MC precision. Moreover similar method was 
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applied to the 2D GL model which was simulated exten- 
sively (Hu and MacDonald, 1993; Hu et ai, 1994; Kato 
and Nagaosa, 1993; O'Neill and Moore, 1993; Tesanovic 
and Xing, 1991) and for which longer series are available 
(Brezin et ai, 1990; Hikami et ai, 1991) and agreement 
is still perfect. We conclude that the method is precise 
enough to study the melting problem. 

Now let us mention several issues, which prevented 
its use and acceptance early on. Ruggeri and Thou- 
less (Ruggeri and Thouless, 1976; Ruggeri, 1978; Thou- 
less, 1975) tried to use BP to calculate the specific heat 
without much success because their series were too short 
(Wilkin and Moore, 1993). In addition they tried to force 
it to conform to the solid expression at low temperatures, 
which is impossible. Attempts to use BP for calculation 
of melting also ran into problems. Hikami, Fujita and 
Larkin (Brezin et ai, 1990; Hikami et ai, 1991) tried to 
find the melting point by comparing the BP energy with 
the one loop solid energy and obtained ot = —7. How- 
ever their one loop solid energy was incorrect (by factor 
v2) and in any case it was not precise enough, since the 
two loop contribution is essential. 

To conclude the BP method and the OPT are pre- 
cise enough to quantitatively determine thermodynamic 
properties of the vortex liquid, including the supercooled 
one. The precision is good enough in order to determine 
the melting line. We therefore turn to the physical conse- 
quences of the analytical methods for both the crystalline 
and the melted liquid states. 

Magnetization and specific heat in vortex liq- 
uids 

As long as the free energy is known, one differentiates 
it to calculated other physical quantities like entropy, 
magnetization and specific heat using general LLL for- 
mulas. Since the BP formulas, although analytical, are 
quite bulky (and can be found in Mathematica file) we 
will not provide them. The magnetization curves were 
compared to those in fully oxidized YBCO in (Li and 
Rosenstein, 2002a) to data of ref. ( Nishizaki et ai, 2000) 
and with Nb (after correction to a rather small n) to data 
of ref. (Salem-Sugui et ai, 2002), while the specific heat 
data were compared with experimental in SnNb3 in ref. 
(Lortz et ai, 2006) and in Nb (also after the finite k 
correction) . 



E. First order melting and metastable states 

1. The melting line and discontinuity at melt 
Location of the melting line 

Comparing solid two - loop free energy given by 
eq.(188) and liquid BP energy, Fig. 11, we find that they 
intersect at a™ = —9.5 (see insert for the difference). The 
available 3D Monte Carlo simulations (Sasik and Stroud, 
1995) unfortunately are not precise enough to provide an 
accurate melting point since the LLL scaling is violated 
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Scaled Temperature 

FIG. 11 The melting point and the spinodal point of the 
crystal. Free energy of the crystalline and the liquid states 
are equal at melt, while metastable crystal becomes unstable 
at spinodal point. 

and one gets values a™ = —14.5,-13.2,-10.9 at magnetic 
fields 1T,2T,5T respectively. This is perhaps due to small 
sample size (~100 vortices). The situation in 2D is bet- 
ter since the sample size is much larger. We performed 
similar calculation to that in 3D for the 2D LLL GL liq- 
uid free energy, combined it with the earlier solid energy 
calculation (Li and Rosenstein, 2002a; Rosenstein, 1999) 

^ = -#+21ogM_l^_ 2 . 92 . (248) 
vol 2(3 A 47r 2 a\ 

and find that the melting point obtained a™ = —13.2. 
It is in good agreement with numerous MC simulations 
(Hu and MacDonald, 1993; Hu et ai, 1994; Kato and 
Nagaosa, 1993; Li and Nattermann , 2003). 

Comparison with phenomenological Lindemann 
criterion and experiments 

Phenomcnologically melting line can be located using 
Lindemann criterion or its more refined version using De- 
bye - Waller factor. The more refined definition is re- 
quired since vortices are not point - like. It was found 
numerically for Yukawa gas (Stevens and Robbins, 1993) 
that the Debye - Waller factor e~ 2W (ratio of the struc- 
ture function at the second Bragg peak at melting to its 
value at T = 0) is about 60%. To one loop order one gets 
using methods of (Li and Rosenstein, 1999b) to calculate 
the Debye - Waller factor at the melting line obtained 
here by using the non-perturbative method 

e~ 2W = 0.50. (249) 

The higher loop correction to this factor is supposed to 
be positive and the total value might be equal to a value 
around 0.6 (we did not undertake this calculation due to 
the complexity). However we apply a "one loop" criterion 
(Debye - Waller factor is 0.5 calculated to one loop), and 
this method was applied to the layered superconductor 
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based on Lawrence- Doniach- Ginzburg Landau model 
and the rotating Bose Einstein condensate in (Feng et al, 
2009; Wu et al, 2007), and the results were both in sur- 
prisingly good agreements with numerical calculations in 
(Cooper et al, 2001; Hu and MacDonald, 1997). 

The melting line in accord with numerous experiments 
in both clean low T c materials like NbSe2( Adesso et al, 
2006; Kokubo et al, 2004, 2005; Kokubo et al, 2007; 
Thakur et al, 2005; Xiao et al, 2004) and Nb 3 Sn(Lortz 
et al, 2006), in which the line can be inferred from the 
peak effect (see below) and various dynamical effects or 
high T c like the fully oxidized YBa 2 Cu 3 7 ( Nishizaki 
et al, 2000), see fit in (Li and Rosenstein, 2002a). The 
fully oxidized YBCO is best suited for the application of 
the present theory, since pinning on the mesoscopic scale 
is negligible. For example the melting line is extended 
beyond 30T as shown in (Li and Rosenstein, 2002a). 

Melting lines of optimally doped untwined (Bouquet 
et al, 2001; Schilling et al, 1996, 1997; Welp et al, 
1991, 1996; Willcmin et al, 1998) YBa 2 Cu 3 7 ^ s and 
DyBa 2 Cu 3 7 (Revaz et al, 1998; Roulin et al, 1996a, b) 
are also fitted extremely well(Li and Rosenstein, 2003). 
More recently both NbSe2 and thick films of Nb 3 Ge were 
fitted in ref. ( Kokubo et al, 2007) in which disorder 
is significant, but the pristine melting line is believed 
to be clearly seen in dynamics via peak effect. In Ta- 
ble 2 parameters inferred from these fits are given where 
the data for YBC0 7 ^ s , YBC0 7 , DyBC0 6 . 7 are taken 
from (Schilling et al, 1996), ( Nishizaki et al, 2000), 
(Roulin et al, 1996a) respectively. Parameters like Gi 
characterizing the strength of thermal fluctuations dif- 
fer a bit from the often mentioned (Blatter et al, 1994). 
Similar fits were made in 2D for organic superconductor 
(Fruchter et al, 1997). Unlike the Lindcmann criterion, 
the quantitative calculation allows determination of vari- 
ous discontinuities across the melting line (since we have 
energies of both phases) to which we turn next. 

Table 2 Parameters of high T c superconductors 
deduced from the melting line 



material 


T c 


H C 2 


Gi 


K 


7a 


YBC0 7 s 


93.07 


167.53 


1.910~ 4 


48.5 


7.76 


YBC0 7 


88.16 


175.9 


7.010~ 5 


50 


4 


DyBC0 6 . 7 


90.14 


163 


3.210~ 5 


33.77 


5.3 



2. Discontinuities at melting 
Magnetization jump 

The scaled magnetization is defined by m (a-r) = 
— ^t/(ot) can be calculated in both phases and the dif- 
ference Am = m s — mi at the melting point a™ = —9.5 
is 



This was compared in ref. (Li and Rosenstein, 2003) with 
experimental results on fully oxidized Y Ba 2 Cu 3 7 ( 
Nishizaki et al, 2000) and optimally doped untwined 
YBa 2 Cu 3 7 s (Welp et al, 1991, 1996). These samples 
probably have the lowest degree of disorder not included 
in calculations. 

Specific heat jump 

In addition to the delta function like spike at melting 
following from the magnetization jump discussed above 
experiment shows also a specific heat jump (Bouquet 
et al, 2001; Lortz et al, 2006; Schilling et al, 1996, 1997). 
The theory allows to quantitatively estimate it. The spe- 
cific heat jump is: 



2 - 2b + t 
t 



Ac = 0.0075 
-0.20G* 1 / 3 (b-l-t) (A) 



(251) 



2/3 



It was compared in ref. (Li and Rosenstein, 2003), 
with the experimental values of ( Willcmin et al, 1998). 
See also comparison with specific heat in NbSn 3 of ref. 
(Lortz et al, 2006). 

In addition the value of the specific heat jump in the 
2D GL model is in good agreement with MC simulations 
(Hu and MacDonald, 1993; Hu et al, 1994; Kato and Na- 
gaosa, 1993), while the 3D MC result is still unavailable. 



3. Gaussian approximation in the crystalline phase and the 
spinodal line 

Gaussian Variational Approach with shift of 
the field 

Gaussian variational approach in the phase exhibiting 
spontaneously broken symmetry is quite a straightfor- 
ward, albeit more cumbersome, extension of the method 
to include " shift" v (r) . In our case of one complex field 
one should consider the most general quadratic form 

K= [ [r{r)-v*{r)]G-\r,r')W) - v(r')} 

J r,r' 

+ [ip(r) - v(r)} H* (r, r') [ip{r') - v{r')\ + c.c (252) 

To obtain "shift" v and "width of the gaussian" which 
is a matrix containing G and H , one minimizes the gaus- 
sian effective free energy (Cornwall et al, 1974), which 
is an upper bound on the energy. Assuming hexagonal 
symmetry (a safe assumption for the present purpose), 
the shift should be proportional to the zero quasi - mo- 
mentum function, v (r) = vip(r), with a constant v taken 
real thanks to the global U(l) gauge symmetry. On LLL, 
as in perturbation theory, we will use the phonon vari- 
ables Ok and Ak defined in quasimomentum basis eqs. 
(139), (143) instead of V(r) 



AM _ Am 

~m7 = 



ip(r) = vip(r) + 



1 



= 0.018 



(250) 



V2(27rp 



ik z z 



Ck^k(r) (O k 



iA k ) . 
(253) 
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The phase defined after eq. (148) is quite important for 
simplification of the problem and was introduced for fu- 
ture convenience. The most general quadratic form in 
these variables is 

K= [ O k Goo(k)0- k + A k G- A \{k)A_ k (254) 

Jk 

+O k G Q \{k)A_ k + A k G OA (k)0_ k , 

with matrix of functions to be determined together with 
the constant v by the variational principle. The gaussian 
free energy is 



gauss 

vol 



2 , Pa 4 
a T v + —v 



2 5 /V 



1 



2 (2nY 



1 



-{(k 2 z /2 + a T ) 



logtdet^)],^ 

x [Goo (k) + Gaa (k)} + v 2 [(2/? k + | 7k |) G QO (k) 

1 



(2p k -\ 7k \)G AA (k)}} 



4(2tt) c 



(255) 



{-U I | 7 k| {Goo (k) Gaa {k))f + 4[ / | 7k | G OA (k)} 2 
2/3a J k h 

Pk-i [Goo (k) + G A a (k)} [G 00 (I) + G AA (I)}}, 



k,l 

leading to the following minimization equations are: 
a,T 1 



v 2 + 



Pa 



2 (2tt) 3 /3 a 



(2/3 k +| 7 k |) Goo (k) 



+ (2/3 k -|7k|) Gaa (fc) 



2 5 /\ [G(k)- 1 ] OQ = k 2 z /2 + a T + v 2 (2/3 k + | 7k |) + 



/K 1 



l7k| Iti 

Pa 



Goo (I) + 2/? k -i - 



l7k| hj 

(3a 



implicitly appears on the left hand side due to a need to 
invert the matrix G. Obviously GoA(k) = is a solution 
and in this case the matrix diagonalizcs. However gen- 
eral solution can be shown to differ from this simple one 
just by a global gauge transformation. Subtracting the 
OO equation from the AA equation above, eq. (256) and 
using the OA equation, we observe that matrix G _1 has 
a form: 



G- 1 



( Gon{k) G~ AO (k) 



J oo 

G- AO {k) G- A \{k) 



1 k 2 /2 + „ 2 0k 



AO^ 
AA ( 
2 /O _L ,,2 



(258) 



2 5 /2tT 



MAOk 



MAOk 

k 2 /2 + M^ k 



with 



2 

Mok 

2 

MAOk 



= £" k + Ai | 7k | ; /4 k = E(k) - Ai | 7 k 
= A 2 |7k| 



(259) 



where Ai,A2 are constants. Substituting this into the 
gaussian energy one find s that it depends on Ai, A 2 via 
the combination A = \/ A 2 + A\ only. Therefore with- 
out loss of generality we can set A 2 = 0, thereby return- 
ing to the Goa = case. 

Using this observation, the gap equations significantly 
simplify. The function £? k and the constant A satisfy: 



E k = a T + 2v 2 (3 k + 2x 



(256) 
1 

2 (2tt) 3 
Gaa (0 



, — + 
VMoi Mai 



< /3 k - 1 
/3a A = -or - 2 < (3i 

and shift equation 



>i; 



-+-) 

MOl MAI / 



(260) 



>i. (261) 



and 

l 5/2 *[G(k)-i] AA 
2/? k _ 1+ l7k|1711 



| + ar + , 2(2/3fe -| 7k |) + -J_, 



(3a 

2V\[G{kTX A 
J7k| 



Gaa (0 + 2/3 k _! 



l7k| Iti 



Pa 

2 5 /2 7rGoA (fc) 



Goo (0 



G 00 {k)G A A{k) - G OA {kf 



Pa 2(2tt) 



|7i| Goa (I) 



(257) 



These equations look quite intractable, however they 
can be simplified. 

How to eliminate the off - diagonal terms 

The crucial obs ervation is that after we have inserted 
the phase c k = ^/ 7k / | 7k | in eq. (255) using our expe- 
rience with perturbation theory, Gao appears explicitly 
only on the right hand side of the last equation. It also 



Pa 



< 



2/3 k 



I7k| 



2/3 k 



1 7k | 



MOk 



MAk 



> 



(262) 



The gaussian energy (after integration over k z ) be- 
comes: 



J_ 

vol 



"V + y/+/!+/ 2 +/ 3 ; 



A =< MOk + MAk >k5 

f 2 =a T < (fi Q l + n A l) + v 2 [(2(3 k + | 7 fc|) M ok 

+ (2(3 k - \ lk \) v A l] > k ; (263) 
h =< Pk-i (MoL + VaI) (Moi + MaO >k,i 
+ 2/t ^ ^^ol- MaO >kf- 

The problem becomes quite manageable numerically af- 
ter one spots an unexpected small parameter. 
The mode expansion 
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Using a formula eq.(416) 

oo 

Pk = ^X"ft(k) (264) 

n=0 

/3„(k) = ^ exp[ik.X], 

|X| 2 =rmi 

derived in Appendix A and the hexagonal symmetry of 
the spectrum, one deduces that Eu can be expanded in 
"modes" 

£ k = ^£„/3„(k) (265) 

The integer n determines the distance of a points on re- 
ciprocal lattice from the origin, and \ = e x p[ — °a/2] = 
cxp[— 27r/v / 3] = 0.0265. One estimates that E n ~ x™°t, 
therefore the coefficients decrease exponentially with n. 
Note that for some integers, for example n = 2,5,6, 
(3 n = 0. Retaining only first s modes will be called "the 
s mode approximation". We minimized numerically the 
gaussian energy by varying v, A and first few modes of 



Table 3. 

Mode expansion. 





-30 


-20 


-10 


-5.5 


/ 


-372.2690 


-159.5392 


-33.9873 


-6.5103 



The sample results of free energy density for various ar 
with 3 modes are given in Table 3. In practice two modes 
are also quite enough. We see that in the interesting 
region of not very low temperatures the energy converges 
extremely fast. In practice two modes are quite enough. 

Spinodal point 

One can show that above 

a sp t nodal = _ 5 5 (2g6) 

there is no solution for the gap equations. The corre- 
sponding value in 2D is a s P modal — _7 anc i j s consistent 
with the relaxation time measured in Monte Carlo simu- 
lations ref.(Kato and Nagaosa, 1993). The spinodal point 
was observed in NbSe 2 ( Adesso et ai, 2006; Thakur 
et ai, 2005; Xiao et ai, 2004) at the position consistent 
with the theoretical estimate. 

Corrections to the gaussian approximation 
The lowest order correction to the gaussian approx- 
imation (that is sometime called the post - gaussian 
correction) was calculated in ref. (Li and Rosenstein, 
2002a, b,c) to determine the precision of the gaussian ap- 
proximation. This is necessary in order to fit experiments 
and compare with low temperature perturbation theory 
and other nonperturbative methods. 

A general idea behind calculating systematic correc- 
tions to the gaussian approximation was already de- 
scribed for liquid in subsection C and modifications are 



quite analogous to those done for the gaussian approxi- 
mation. Results for the specific heat were compared in 
ref. (Li and Rosenstein, 2002c). Generally the post - gaus- 
sian result is valid till ar = — 7 and rules out earlier ap- 
proximations, as the one in ref . (Tesanovic et ai, 1992; 
Tesanovic and Andreev, 1994) (dotted line). 



IV. QUENCHED DISORDER AND THE VORTEX GLASS. 

In any superconductor there are impurities cither 
present naturally or systematically produced using the 
proton or electron irradiation. The inhomogeneities both 
on the microscopic and the mesoscopic scale greatly af- 
fect thermodynamic and especially dynamic properties 
of type II superconductors in magnetic field. The field 
penetrates the sample in a form of Abrikosov vortices, 
which can be pinned by disorder. In addition, in high T c 
superconductors, thermal fluctuations also greatly influ- 
ence the vortex matter, for example in some cases ther- 
mal fluctuations will effectively reduce the effects of dis- 
order. As a result the T — H phase diagram of the high 
T c superconductors is very complex due to the competi- 
tion between thermal fluctuations and disorder, and it is 
still far from being reliably determined, even in the best 
studied superconductor, the optimally doped YBCO su- 
perconductor. 

It is the purpose of this section to describe the glass 
transition and static and thermodynamic properties of 
both the disordered reversible and the irreversible glassy 
phase. The disorder is represented by the random com- 
ponent of the coefficients of the GL free energy, eq.(20), 
and the main technique used is the replica formalism. 
The most general so called hierarchical homogeneous (liq- 
uid) Ansatz (Mezard, 1991) and its stability is considered 
to obtain the glass transition line and to determine the 
nature of the transition for various values of the disor- 
der strength of the GL coefficients. In most cases the 
glassy phase exhibits the phenomenon of "replica sym- 
metry breaking, when ergodicity is lost due to trapping 
of the system in multiple metastable states. In this case 
physical quantities do not possess a unique value, but 
rather have a distribution. We start with the case of 
negligible thermal fluctuations. 



A. Quenched disorder as a perturbation of the vortex 
lattice 

1. The free energy density in the presence of pinning potential 
GL model with ST C disorder 

We start with space variations of the coefficient of | \I/ 1 2 , 
eq.(20) distributed as a white noise, eq.(21). It can be 
regarded as a local variation of T c . As was mentioned in 
section I other types of disorder are present and might 
be important, however, as will be shown later are more 
complicated. 
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Since a point - like disorder breaks the translational 
symmetry in all directions including that of the magnetic 
field z, one has to consider configurations dependent on 
all three coordinates and take into account anisotropy, 
discussed in subsection IE. We restrict to the case m* — 



F[W] = 



2m* 



iD*r + 



2m* 



\d z t>\ 2 + a(T-T c 



[l + W{r)\ |1f 



(267) 



where W(r) is the 5T C random disorder (real) field, which 
we assume to be a white noise with variance that can be 
written in the following form: 



W(r)W(r>) = nfCcS 3 (r - r') . 



(268) 



The dimensionless parameter n is proportional to the 
density of pinning centers and a single pin's "strength", 
while £ c = £ (m* /m*) 1 ^ 2 is the coherence length in the 
field direction. The units we use here are the same as 
before with the addition of £ c as the unit of length in 
the z direction. As in previous sections, we will con- 
fined ourselves mainly to the region in parameter space 
described well by the lowest Landau level approximation 
(LLL) defined next. 

The disordered LLL GL free energy in the 
quasi-momentum basis 

In the units and the field normalization described in 
IIA the LLL energy becomes: 



F[W]= J[\\d^\ 2 -a H \^\ 
+ i_^ r) |vl,| 2 + i|vl/| 4 ], 



(269) 



where cm = h (1 — b — t) and 



W(r)W(r') = n5 6 (r - r') 



(270) 



in the new length unit. The order parameter field on LLL 
can be expanded in the quasi - momentum basis defined 
in IIIA as 



(271) 



where k = (k, fc z ), functions are defined in 
eqs.(134),(137) and the integration measure was defined 
in section IIIA to be the Brillouin zone in the x — y plane 
and the full range of momenta in the z direction. We 
consider the hexagonal lattice, although modifications re- 
quired to consider a different lattice symmetry are minor. 
Using the quasimomentum LLL functions of eq.(134), the 
disorder term becomes 



F di . 



l-t 



W(r)\V(r)\ z = f n.,vT/;tt 

Jk.i 



(272) 



with 



(273) 



(274) 



l-t f 

Wk < 1 = TTTTTs / W M M W (r) ■ 
The rest of the terms can be written as 

F clean - f {k 2 z /2 - a H ) Vy k + — x 

[ [k, k'\i, i'] *y* k , $,fTi(H^' -i -i' - Q) 

Jk,k',l,V q 

with [k,l\k'l']^= ^j;^(r) W (r)^(r) W (r) and 

where Q = , 0^ and Q is the reciprocal lattice vectors 

as k,l,k',l' satisfy the momentum conservation up to a 
reciprocal lattice vector, [fc, l\k'l'] will be equal to zero if 
k + k' -l-l' ^Q. 



2. Perturbative expansion in disorder strength. 

Expansion around the Abrikosov solution 

The GL equations derived from the free energy in the 
quasimomentum basis are 

(k 2 z /2 -a H )V k +a ( w k . l ^ l + [ (275) 

Jl Jk'l,V 

J2S(k + k'-l-l'-Q) % *„ = 0. 



Q 



(2nT 



The parameter a = 1 inserted there will help with count- 
ing orders. The expansion in orders of the disorder 
strength a reads: 



(276) 



The clean case Abrikosov solution of section II is defined 
as the quasimomentum zero. Therefore 



= (2tt) 



3/2 




5k- 



(277) 



The delta function appears due to its long - range trans- 
lational order. Now the equation eq.(275) can be solved 
order by order in a. Since contributions linear in dis- 
order potential will average to zero, in order to get the 
leading contribution of disorder one should calculate the 
free energy to the second order in a. Multiplying exact 
equation eq.(275) by \&* and integrating over k, one can 
express the order four in * term via simpler quadratic 
ones: 



k 2 z /2-a H )\* k 



+ - / (278) 
z Jk,i 
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Substituting the expansion eq.(276) and using delta func- 
tions of *(°) of eq.(277) one gets the following a 2 terms 



J?(2) 



3/2 ,„ n3/2 

a H ( 2?r ) r,T,( 2 ) 



$(1) 



2/?i /2 

2 + 4 /2 (27T) 



2/3i /2 



[*r+^ 2) ] + i/(^/2-^) 

/K fc *i 1) + *i 1 >*«; M ]. 
Jk 



(279) 



Therefore the second order correction to is needed only 

for zero quasi - momentum. 

First order elastic response of the vortex lattice 
To order a one obtains the following equation 



and using the expression for eq.(281) one obtains 
various terms quadratic in disorder w. The disorder av- 
erages are 



Wk,lWk',l' = 



(1-tfnV 



[k,k'\l,l'}; 



w k,iw k ',i 
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[k,l'\l,k'}; (285) 



and so that the pinning energy becomes after some alge- 
bra 



(fc 2 /2 - a H ) + w kfi (2*f 2 J%L+ (280) 

V Pa 

2^A*i 1) + ^7 k * ( ? = 

PA PA 

as Q = because of the conservation of quasimomentum 
in this case. This equation and its complex conjugate 
lead to a system of two linear equations for two vari- 
ables \E , ( 1 ) and Solution, not surprisingly, involves 
the spectrum of harmonic excitations of the vortex lattice 
already familiar from the perturbative corrections due to 
thermal fluctuations, IIIA: 



, x3/2 1/2 

£ k £ k P A ' PA 

a H * i 
—o-1kW_ kQ \, 
PA 

where e£, e° are defined in eq.(150). 

Disorder average of the pinning energy to lead- 
ing order 

The relevant equation (zero quasi - momentum) at the 
second order in a is: 

-a H ^ + [ m , k *W + a -f \20 A *W 

Jk PA 1 J 

1/2 „ 

(3 A /2 (2nf 2 Jl 
leading to 

"*< 2 > + *( 2 > =- / *< 1 >{«;o,fc+ (288) 
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Substituting this into eq.(279) and simplify the equa- 
tion using eq.(280), we obtain the energy expressed via 

k 

F {2) = ""Ifr I + ^.o} (284) 
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Integrating over k z , one obtains finally, 



(l-t) 2 na H < P h +\ 2 J + ^hJ > (2g7) 
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where are given in eq.(150). 

Using expansion for small k of the functions /3k and 
|7 k | derived in Appendix A, one can see that the second 
term is finite 



«^ 2 k. 



Numerically 
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Stronger disorder: 2D GL and columnar de- 
fects 

The same calculation can be performed in 2D with the 
result 



H-tYna H / ( A^ + A i _b4 (290) 



vol 
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This is logarithmically IR divergent at any value of the 
disorder strength 



d 2 k. 



(291) 



Therefore either the dependence is not analytic or (more 
probably) disorder significantly modifies the structure of 
the solution. Generalization in another direction, that 
of long - rage correlated disorder can also be easily per- 
formed. One just replaces the white noise variance by a 
general one 



W(r)W(r') =K(r-r'). 



(292) 



38 



For columnar defects the variance is independent of z, 
K(r-r')=nS(r-r'), (293) 
and one again obtains a logarithmic divergence. 



vX.N, rj | i-) ^ > 



vol 
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3. Disorder influence on the vortex liquid and crystal. Shift of 
the melting line 

Disorder correction to free energy 

Thermal fluctuations in the presence of quenched are 
still described by partition function 



1 



F[*] 



1 - 1 
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If W is small, we can calculate Z by perturbation theory 
in W. To the second order free energy — TlnZ is 



G = G 
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where () and f dean denote the thermal average and free 
energy of the clean system. Averaging now over disorder 
one obtains 



AGrf,-., — G — G 
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Therefore one has to calculate the superfluid density ther- 
mal correlator. In LLL approximation and LLL units, 



AG dls /V = -r(t)Af; 
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r(i) = -^(l-<) 2 = 
4w t 
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Calculations of this kind in both solid and liquid were 
the subject of the previous section. 

Correlators in the crystalline and the liquid 
states 

Within LLL (and using the LLL units introduced in 
section III) the one loop disorder correction to the crys- 
tal's energy is: 



crystal 



2.U\a T 



1/2 
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An explicit expression for fn q (aT), obtained using the 
Borel-Pade resummation of the renormalized high tem- 
perature series (confirmed by optimized Gaussian series 




FIG. 12 Phase diagram for YBCO. 



and Monte Carlo simulation) is rather bulky and can be 
found in ref. (Li and Rosenstein, 2002a). One can de- 
rived an expression for the disorder correction in liquid 
by differentiating the "clean" partition function with re- 
spect to parameters: 



Hi 
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These two results enable us to find the location of the 
transition line and, in addition, to calculate discontinu- 
ities of various physical quantities across the transition 
line. 

The "downward shift" of the first order transi- 
tion line in the T — H plane 

It was noted in section III that in a clean system a 
homogeneous state exists as a metastable overcooled liq- 
uid state all the way down to zero temperature (not 
just below the melting temperature corresponding to 
ax = —9.5, see the n = line in Fig. 12) . This is of 
importance since interaction with disorder can convert 
the metastable state into a stable one. Indeed generally 
a homogeneous state gains more than a crystalline state 
from pinning, since it can easier adjust itself to the topog- 
raphy of the pinning centers. At large |ar| in particular 
Af liq oc a\ compared to just A/ so/ cx IotI 1 ^ 2 - As a result 
in the presence of disorder the transition line shifts to 
lower fields. The equation for the melting line is 

d{a T ) = (fu q ~ f S oi)/(Afu q - A/ so/ ) = n(t). (301) 

The universal function d{ax), plotted in Fig. 13 turns out 
to be non-monotonic. This is an important fact. Since 
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FIG. 13 Universal function d{ar) determining the shift of 
the melting line due to disorder 



n(t) is a monotonic function of t, one obtains the tran- 
sition lines for various n in Fig. 12 by "sweeping" the 
Fig. 13. A peculiar feature of d{ax) is that it has a lo- 
cal minimum at ax ~ —17.2 and a local maximum at 
ax ~ —12.1 (before crossing zero at ax ~ — 9.5). There- 
fore between these two points there are three solutions to 
the melting line equation. As a result, starting from the 
zero field at T c , the transition held H{T) reaches a maxi- 
mum at E ent beyond which the curve sharply turns down 
(this feature was called "inverse melting" in (Avraham 
et al., 2001)) and at E mag backwards. Then it reaches a 
minimum and continues as the Bragg glass - vortex glass 
line roughly parallel to the T axis. 

The temperature dependence of the disorder strength 
n(t), as of any parameter in the GL approach should 
be derived from a microscopic theory or fitted to experi- 
ment. General dependence near T c is: n(t) = n(l~t) 2 /t . 
The extra factor (1 — t) 2 , not appearing in a phenomcno- 
logical derivation (Blatter et ai, 1994), is due to the fact 
that near H c i order parameter is small oc (1 — t) and 
disorder (oxygen deficiencies) locally destroys supercon- 
ductivity rather than perturbatively modifies the order 
parameter. The curves in Fig. 12 correspond to the dis- 
order strength n = 0.08,0.12,0.3. The best fit for the 
low field part of the experimental melting line H m (T) of 
the optimally doped YBCO (data taken from (Schilling 
et ai, 1996), T c = 92.6, 7 = 8.3) gives Gi = 2.0 x 10~ 4 , 
H C 2 = 190 T, k = A/£ = 50 (it is consistent with other 
experiments, for example, (Dcligiannis et al, 2000; Shi- 
bata et ai, 2002)). This part is essentially independent 
of disorder. The upper part of the melting curve is very 
sensitive to disorder: both the length of the "finger" and 
its slope depend on n . The best fit is no = 0.12. This 
value is of the same order of magnitude as the one ob- 
tained phcnomcnologically using eq.(3.82) in ref. (Blat- 
ter et ai, 1994). We speculate that the low temperature 
part of the "unified" line corresponds to the solid - vor- 
tex glass transition H* (T) observed in numerous experi- 
ments (Bouquet et ai, 2001; Kokkaliaris et ai, 2000; Pal 
et al, 2001, 2002; Radzyner et ai, 2002; Schilling et al, 



1996, 1997; Shibata et ai, 2002), see data (squares in 
Fig. 12) taken from ( Shibata et ai, 2002). A complicated 
shape of the "wiggling" line has been recently observed 
(Pal et ai, 2001, 2002). Now we turn to a more detailed 
characteristics of the phase transition. 

Discontinuities across the transition and the 
Kauzmann point. Absence of a second order 
transition. 

Magnetization and specific heat of both solid and liq- 
uid can be calculated from the above expressions for free 
energy. Magnetization of liquid along the melting line 
H m {T) is larger than that of solid. The magnetization 
jump is compared in (Li and Rosenstein, 2003)a with the 
SQUID experiments (Schilling et al., 1997) in the range 
80 — 90K (triangles) and of the torque experiments (stars 
( Willemin et ai, 1998) and circles ( Shibata et ai, 2002) 
). One observes that the results of the torque experi- 
ments compare surprisingly well above 83K , but those of 
( Shibata et ai, 2002) vanish abruptly below 83K unlike 
the theory and are inconsistent with the specific heat ex- 
periments (Deligiannis et ai, 2000; Schilling et ai, 1996) 
discussed below. The SQUID data are lower than the- 
oretical (same order of magnitude though). We predict 
that at lower temperatures (somewhat beyond the range 
investigated experimentally so far) magnetization reaches 
its maximum and changes sign at the point E mag (at 
which magnetization of liquid and solid are equal) . 

In ref. (Li and Rosenstein, 2003), entropy jump was 
calculated using the Clausius - Clapeyron relation 
dH m {T)/dT = — AS/ AM and compared with an exper- 
imental one deduced from the spike of the specific heat 
((Schilling et al., 1996), and an indirect measurement 
from the magnetization jumps (in ref.( Shibata et al, 
2002)). At high temperatures the theoretical values are 
a bit lower than the experimental and both seem to ap- 
proach a constant at T c . The theoretical entropy jump 
and the experimental one of (Schilling et ai, 1996) vanish 
at E ent (Fig. 12) near 75K. Such points are called Kauz- 
mann points. Below this temperature entropy of the liq- 
uid becomes smaller than that of the solid. Note that the 
equal magnetization point E mag is located at a slightly 
lower field than the equal entropy point E ent . Experi- 
mentally a Kauzmann point was established in BSCCO 
as a point at which the "inverse melting" appears (Avra- 
ham et al., 2001). The Kauzmann point observed at a 
lower temperature in YBCO in ref. (Radzyner et al, 
2002) is different from E ent since it is a minimum rather 
than a maximum of magnetic field. It is also located 
slightly outside the region of applicability of our solu- 
tion. The point E ent is observed in (Pal et ai, 2001, 
2002) in which the universal line is continuous. 

In addition to the spike, the specific heat jump has also 
been observed along the melting line H m (T) (Deligiannis 
et ai, 2000; Schilling et ai, 1996, 1997). Theoretically 
the jump does not vanish either at E ent or E mag , but 
is rather flat in a wide temperature range. Our results 
are larger than experimental jumps of (Schilling et al., 
1996) (which are also rather insensitive to temperature) 
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by a factor of 1.4 to 2 (Li and Rosenstein, 2003). In 
many experimental papers there appears a segment of the 
second order phase transition continuing the first order 
melting line beyond a certain point. In (Bouquet et ai, 
2001) it was shown that at that point the specific heat 
profile shows "rounding" . We calculated the specific heat 
profile above the universal first order transition line. It 
exhibits a " rounding" feature similar to that displayed by 
the data of (Bouquet et ai, 2001; Schilling et al., 1996, 
1997, 2002) with no sign of the criticality. The height 
of the peak is roughly of the size of the specific heat 
jump. We therefore propose not to interpret this feature 
as an evidence for a second order transition above the 
first order line. 

Limitations of the perturbative approach 

Of course the perturbative approach is limited to small 
couplings only. In fact, when the correction is compared 
to the main part of the lattice energy the range become 
too narrow for practical applications in low temperature 
superconductors. For high T c superconductors thermal 
fluctuations cannot be neglected at higher temperatures 
since it "melts" the lattice and even at low temperatures 
provides thermal dcpinning. On the conceptual side, it is 
clear that disorder contributes to destruction of the trans- 
lational and rotational order. Therefore at certain disor- 
der strength, vortex matter might restore the translation 
and rotation symmetries, even without help of thermal 
fluctuations. It is possible to use the perturbation the- 
ory in disorder with the liquid state as a starting point in 
the case of large thermal fluctuations, however it fails to 
describe the most interesting phenomenon of the vortex 
glass introduced by Fisher (Fisher, 1989; Fisher et ai, 
1991). Therefore one should try to develop non - pertur- 
bative methods to describe disorder This is the subject 
the following sections. 



1. Replica approach to disorder 
The replica trick 

The replica method is widely used to study disordered 
electrons in metals and semiconductors, spin glasses and 
other areas of condensed matter physics and far beyond 
it (Itzykson and Drouffe, 1991). It was applied to vor- 
tex matter in the elastic medium approximation (Bogner 
et ai, 2001; Giamarchi and Le Doussal, 1994, 1995a,b, 
1996, 1997; Korshunov, 1993; Nattcrmann, 1990). In the 
following we describe the method in some detail. 

The main problem in calculation of disorder averages 
is that one typically has to take the average of non - 
polynomial functions of the statistical sum eq.(18): 

Z = J V*V**e-M F W +1 ^fr w W*(r)\ 2 }, (302) 

Most interesting physical quantities are calculated by 
taking derivatives of the free energy which is a loga- 
rithm of Z. Applying a simple mathematical identity 
to represent the logarithm as a small power, log (z) — 
lim„_>o - (z™ — 1), the average over the free energy is writ- 
ten as: 



T= -cj t \im-(Z n - 1). (303) 

The quantity Z n can be looked upon as a statistical sum 
over n identical "replica" fields $> a , a = 1, ...,n: 

b J*>> 

(304) 

where -F[\& a ] is the free energy (in physical units mean- 
time) without disorder. Note that the disorder potential 
enters in the exponent. The disorder measure, consistent 
with variance in eq.(268) is a gaussian. Therefore disor- 
der average is a gaussian integral which can be readily 
performed: 



B. The vortex glass 

When thermal fluctuations are significant the efficiency 
of imperfections to pin the vortex matter is generally di- 
minished. This phenomenon is known as "thermal de- 
pinning" . In addition, as we have learned in section III, 
the vortex lattice becomes softer and eventually melts via 
first order transition into the vortex liquid. The inter - 
dependence of pinning, interactions and thermal fluctua- 
tions is very complex and one needs an effective nonper- 
turbative method to evaluate the disorder averages. Such 
a method, using the replica trick was developed initially 
in the theory of spin glasses. It is more difficult to apply 
it in a crystalline phase, so we start from a simpler ho- 
mogeneous phase (the homogeneity might be achieved by 
both the thermal fluctuations and disorder) and return 
to the crystalline phase in the following subsection. 



— [ VWe,-^^ w2 ^Z n [W] (305) 
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where 



F n = £ F[*„] + \r (t) W |* a | 2 l^l 2 • (306) 



After the disorder average different replicas are no longer 
independent. In LLL limit and units, 



Z n = [ e i^f F ", (307) 

n „ J, J T 



a.b 
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This statistical physics model is a type of scalar field 
theory and the simplest nonpcrturbativc scheme com- 
monly used to treat such a model is gaussian approxima- 
tion already introduced in IIIB. Its validity and precision 
can be checked only by calculating corrections. 

Correlators and distributions 

Correlators averaged over both the thermal fluctua- 
tions and disorder can be generated by the usual trick 
of introducing an external "source" into statistical sum 
eq.(302): 

Z[W,S*,S]= f exp{--[F[Sf] + hl-t) (308) 
x f W{r)\^{r)\ 2 + ( * (r) S* (r) + ** (r) S (r)]} 

J r J r 

e iF[*,ff(r),S(r),S'(r^ 



/ 



and taking functional derivatives of the free energy in the 
presence of sources 

F[W,S*,S\ = -w t \og Z[W,S*,S\. (309) 

The first two thermal correlators are 

(* ( r )) = - f * (r) e -^m,W(r),S(r),S'(r)] 

~Z[W,0,0] ss* {r) Z ^ S ^ S ^=o 

= jJ-rrJ r [W,S*,S}\s,s*= ; (310) 

<** W * (r'))c = (** W * (r 7 )) - (** (r)) (* (/)) 
S 2 

= ss(r)ss* ^rf^ s *^\^=o. 

Now the disorder averages of these quantities are made 
using the replica trick 
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(nr)) = -^j^^-(ZlS,ST-l) (311) 
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a a a 

Similar calculation for the two field correlator result in 
= i E <** a (r) * b (/)> . (312) 



a, 6 



In disorder physics it is of interest to know the disor- 
der distribution of physical quantities like magnetization 
(which within LLL is closely related to the correlator, 
see IIIB). The simplest example is the second moment of 
the order parameter distribution (\&* (r)) (<J> (r')). This 



is harder to evaluate due to two thermal averages. One 
still uses eq.(310) twice: 



(r)) (* (/)) = ^ $1 (r) * 2 (/) 

xe-^I'iTOI-^IftW (313) 

W(r)]- 1 i-F[<S 2 ,W(r)] v 



= lim / *i(r)* 2 (r')e _ ^ F[ * l! 
Z™- 2 [W] = lim / *i (r) * 2 (/) e"^ 



E F[*i,W(r)] 
»=i 



The disorder average leads to (Mezard ei a/., 1987): 

(** (r)> (r')> = lim / *i (r) * 2 (r 7 ) e~^ F " 

= hm / *a(r)* 6 (r')e-^ F — Q a . h (314) 

In case of replica symmetry breaking, the formula 
above shall be written as 



<**(r))<*(rO>= hm —L— ^Q a , fc . (315) 
n— *0,a=£b n — 1 ) ^— ' 

Therefor 



a^b 



(**(r)*(r / )) = (** (r) * (r')) c + <** (r)) <* (r')) 

a 

Disordered LLL theory 

Restricting the order parameter to LLL (eq.(112)) 
by expanding it in quasi - momentum LLL functions 
eq.(271), one obtains the disordered LLL theory. Let 
us also rewrite the model in the same units we have used 
in section III. The resulting Boltzmann factor is 2^)2 n .f 

/ = E tyfc* ( fc '/2 + or) r k + fint m} + his, (317) 

a 

with the disorder term 



fdis — ^ ^ 
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(318) 



a.b 2 ( 2?r ) Jk,l,k',l 

x[k,k'|i,ivrv-r4*^, 

in which [k, k'|l, 1'] was defined in eq.(142). 
2. Gaussian approximation 



Gaussian energy in homogeneous (amorphous) 
phase 
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One can recover the perturbative results of the previ- 
ous subsection and even generalize them to finite temper- 
atures, by expanding in r, however the replica method's 
advantage is more profound when nonperturbative ef- 
fects arc involved. We now apply the gaussian approx- 
imation, which has been already used in vortex physics 
in the framework of the elastic medium approach, (Gia- 
marchi and Le Doussal, 1994, 1995a,b, 1996, 1997; Kor- 
shunov, 1990, 1993) following its use in polymer and dis- 
ordered magnets' physics (Mezard, 1991). As usual, ho- 
mogeneous phases are simpler than the crystalline phase 
considered in the previous subsection, so we start from 
the case in which both the translational and the U (1) 
symmetries are respected by the variational correlator: 



{n^l) = G ab {k z ) = ij^u. 



(319) 



Since the gaussian approximation in the vortex liquid 
within the GL approach was described in detail in section 
III, we just have to generalize various expressions to the 
case of n replicas. The gaussian effective free energy is 
expressed via variational parameter (Li and Rosenstcin, 
2002a, b,c; Mezard, 1991) u ab which in the present case is 
a matrix in the replica space. The bubble and the trace 
log integrals appearing in the free energy are very simple: 
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2 5 / 2 7T 



Ub = 2 [u ^ = 2u ab 



2 5 / 2 7T f 

—3 / [LogG-\k z )] aa = 4u aa + const, (320) 
(2ir) Jk 



where the "inverse mass" matrix u ab was defined. As a 
result the gaussian effective free energy density can be 
written in a form: 

n /„ = £ / [LogG- 1 (k z )]aa+ (321) 

-i-3 / [(fc 2 /2 + or) G(k z ) - I] aa + 4 (u aa f -2rJ2 Kb 

(Z7TJ Jk 

= 2^{^ aa + a T u aa + 2(u QQ ) 2 } -2r^\u ab \ 
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where we discarded an (ultraviolet divergent) constant 
and higher order in n, and for simplicity, r (t) is denoted 
by r. 

Minimization equations 

It is convenient to introduce a real (not necessarily 
symmetric) matrix Q a b, which is in one - to - one lin- 
ear correspondence with Hermitian (generally complex) 
matrix u a b via 



Q ab = Rc[M afc ] + lm[u ab ] 



(322) 



Unlike u ab , all the matrix elements of Q a b are indepen- 
dent. In terms of this matrix the free energy can be 
written as 



f !„ = E ( U ^)aa + ^Qaa + 2 (Q aa f -r^Qab- 



a.b 



(323) 

Taking derivative with respect to independent variables 
Qab, gives the saddle point equation for this matrix ele- 
ment: 



drSab + ^QaaSab - 2rQ a b = 0. 



(324) 



Since the electric charge (or the superconducting phase) 
U{1) symmetry is assumed, we consider only solutions 
with real u ab . In this case u a b = Qab is a symmetric real 
matrix. 

The replica symmetric matrices Ansatz and the 
Edwards - Anderson order parameter 

Experience with very similar models in the theory of 
disordered magnets indicates that solutions of these min- 
imization equations are most likely to belong to the class 
of hierarchical matrices, which will be described in the 
next subsection, We limit ourselves here to most obvious 
of those, namely to matrices which respect the Z n replica 
permutation symmetry 



!p{a)p{b) 



(325) 



for any of n! permutations a — > p(a) . If we include also 
disorder in \ip (r)| 4 term, one will find in the low tempera- 
ture region, replica symmetry is spontaneously breaking 
as soon as the Edwards -Anderson (EA) order param- 
eter is non zero. However we will limit our discussion 
to replica symmetric solution (not considering disorder 
in \ip (r)| 4 ) and think that the glass transition appears 
when the EA order parameter is non zero. This transi- 
tion line from zero EA to non zero EA obtained in the 
following is very near to the replica symmetry breaking 
transition line considering weak disorder in |V'( r )| (Li 
et al, 2006b). We also believe that even without disor- 
der in \ip (r)| 4 term, the replica symmetry is breaking if 
we can solve the model non perturbatively. 

The most general matrix of the replica symmetric so- 
lution has a form 



Qab = Uab = S ab U + (1 - Sab) A 



(326) 



The off diagonal elements are equal to the Edwards - 
Anderson (EA) order parameter A. A nonzero value for 
this order parameter signals that the annealed and the 
que nched averages are generally different. Let us calcu- 
late (r)) (tp (r)) starting from the eq.(314). Using the 
eq.(319) , one obtains in gaussian approximation 



(r (r))^(r)) r = 2A 
where () r contains also space average. 



(327) 
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One can visualize this phase as a phase with locally 
broken U (1) symmetry with various directions of the 
phase at different locations with zero average (ip (r)) =0 
but a distribution of non zero value of characteristic 
width A. Distribution of more complicated quantities will 
be discussed in the last subsection. Here we will refer to 
this state as glass, although in the subsection E it will be 
referred to as the "ergodic pinned liquid" (EPL) distin- 
guished from the "nonergodic pinned liquid" (NPL) in 
which, in addition, the crgodicity is broken. Broken er- 
godicity is related to "replica symmetry breaking", how- 
ever, as we show there, in the present model of the 5T C 
disorder and within gaussian approximation RSB does 
not occur. If the EA order parameter is zero, disorder 
doesn't have a profound effect on the properties of the 
vortex matter. We refer to this state just as "liquid". 

The dynamic properties of such phase are generally 
quite different from those of the non -glassy A (zero EA 
order parameter) phase. In particular it is expected to 
exhibit infinite conductivity (Dorsey et al., 1992; Fisher, 
1989; Fisher et al., 1991). However if u a b is replica sym- 
metric, pinning docs not results in the multitude of time 
scales. Certain time scale sensitive phenomena like var- 
ious memory effects (Paltiel et al, 2000a, b; Xiao et al, 
2002) and the responses to "shaking" (Beidenkopf et al, 
2005) are expected to be different from the case when u a b 
breaks the replica permutation symmetry. 

We show in the following subsection that within the 
gaussian approximation and the limited disorder model 
that we consider (the 5T C inhomogeneity only) RSB does 
not occur. After that is shown, we can consider the re- 
maining problem without using the machinery of hierar- 
chical matrices. 

Properties of the replica symmetric matrices 

It is easy to work with the RS matrices like u a t in 
eq.(326). It has two eigenvalues. A replica symmetric 
eigenvector 
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; A' = u+(n- 1)A ~ u-X, (328) 
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where sub leading terms at small n, were omitted in the 
last line, and the rest of the space (replica asymmetric 
vectors) which is n — 1 times degenerate. For example 
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A ~ u — A, 



(329) 



The counting seems strange, but mathematically can 
be defined and works. Numerous attempts to discredit 
replica calculations on these grounds were proven base- 
less. Note that the two eigenvalues differ by order n terms 



only. Projectors on these spaces are 

/ 1 1 .. 1 \ 
1 1 .. 1 



Ps - 



1 



n 



P A =1-P S ; (330) 



1 1 .. 1 
A'P S + AP A . 



Here X is the unit matrix l) ah . It is easy to invert RS 
matrices and multiply them using this form. For example 



■ 1 = A!- 1 P s + &- 1 Pa. 



(331) 



3. The glass transition between the two replica symmetric 
solutions 

The unpinned liquid and the "ergodic glass" 
replica symmetric solutions of the minimization 
equations 

The minimization equation cq.(324) for RS matrices 
takes a form 

- A'- 2 P S - A~ 2 P A + {a T + 4?7) J (332) 
-2r(A / P s + AP A ) = 

Expressing it via independent matrices I and P s one ob- 
tains: 

[A" 2 - A'" 2 - 2r (A' - A)] P s + (333) 
(- A" 2 + a T + Au - 2r A) 1 = 0. 

To leading order in n (first) the P s equation is 

A (A" 3 - r) = 0. (334) 

This means that there exists a RS symmetric solution 
A = 0. In addition there is a non diagonal one. It turns 
out that there is a third order transition between them. 
The second equation, 



a T + 4u - 2rA = 0, 



(335) 



for the diagonal ("liquid") solution, A; = ui, is just a 
cubic equation: 

- u^ 2 + a T + Aui - 2rui = 0. (336) 

For the non - diagonal solution the first equation eq. (334) 
gives A = r 1 / 3 , which, when plugged into the first equa- 
tion, gives: 



u, 



= i(3r 2 / 3 - flT );A=i(3r 2 / 3 -a T )-r- 1 



The matrix u therefore is 

u 9 ab = r-^5 ab + X. 



/3_ 

(337) 
(338) 
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The two solutions coincide when X g = leading to the 
glass line equation 



-1/3 



(3r - 4) . 



(339) 



Free energy and its derivatives. The third order 
glass line. 

Let us now calculate energies of the solutions. Energy 
of such a RS matrix is given, using eq.(331), by 

lf gauss =Y.^~ lp s + ^ l PA) aa + a T u (340) 

~ n {2A" 1 - A~ 2 A + a T u + 2u 2 - r (A 2 + 2AA) } , 

where leading order in small n was kept. The RS energy 
is 



/; = + 2a T ui + 2 (4 - 2r) uf , 



(341) 



which can be further simplified by using eq.(336), 

f l =Aui 1 . (342) 



The glass free energy is even simpler: 

4 = 6^-1 ^ r 2/3_ a ^ 2 



(343) 



Since in addition to the energy, the first derivative of 
the scaled energy, the scaled entropy 



dciT 



= 2r" 1 /s 



and the second derivative, specific heat, 



datp 



(344) 



(345) 



respectively, coincide for both solution on the transition 
line defined by eq.(339). The third derivatives are differ- 
ent so that the transition is a third order one. 

Hessian and the stability domain of a solution 
Up to now we have found two homogeneous solutions of 
the minimization equations. There might be more and 
the solutions might not be stable, when considered on 
the wider set of gaussian states. In order to prove that 
a solution is stable beyond the set of replica symmetric 
matrices u, one has to calculate the second derivative of 
free energy (so called Hessian) with respect to arbitrary 
real matrix Q ab defined in eq.(322): 



n S 2 f eff 



(346) 



{ab)(cd) ~ 2 SQ ab SQ cd 

= i\ [(«- 2 L ( u_1 ) d 6 - i (-- 2 ) ad K 1 ) J + ^(O™ x 

( U ~ 2 )db ~ i ( U_1 )ad ( U ~ 2 )cJ + C - C -} + ^acSbdSab ~ 2r6 ac S bd 



The Hessian should be considered as a matrix in a space, 
which itself is a space of matrices, so that Hessian's in- 
dex contains two pairs of indices of u. We will use a sim- 
plified notation for the product of the Kronecker delta 
functions with more than two indices: 5 ac Sbd5ab = $abcd- 
It is not trivial to define what is meant by "positive def- 
inite" when the number of components approaches zero. 
It turns out that the correct definition consists in finding 
all the eigenvalues of the Hessian "super" matrix. 

Stability of the liquid solution 

For the diagonal solution the Hessian is a very simple 
operator on the space of real symmetric matrices: 



H, 



(ab)(cd) 



Cllabcd + CjJabcd, 



(347) 



where the operators I (the identity in this space) and J 
are defined as 

I = S ac 5 bd ; J = S abcd (348) 
and their coefficients in the liquid phase are: 



cj = 2 (t 



(349) 



with ui being a solution of eq.(336). The corresponding 
eigenvectors in the space of symmetric matrices are 



v^) = A5 cd + B. 



(350) 



To find eigenvalues h of H we apply the Hessian on a 
vector v. The result is (dropping terms vanishing in the 
limit n — > 0): 

H(ab)(cd)Vcd = A (Cj + Cj) S ab + B (c 7 + Cj5 ab ) 

= h(AS ab + B) (351) 

There two eigenvalues are therefore: M 1 ) = cj and = 
cj + cj. Since cj = 4 > 0, the sufficient condition for 
stability is: 



d = 2 (u~ 3 - r) > 0. 



(352) 



It is satisfied everywhere below the transition line of 
eq.(339). 

The stability of the glass solution 

The analysis of stability of the non - diagonal solution 
is slightly more complicated. The Hessian for the non - 
diagonal solution is: 



H, 



(ab)(cd) 



c v V + c v U + cjJ, 



(353) 



where new operators are 

V(ab){cd) = &ac + Sbd', ^(ab)(cd) = 1 (354) 

and coefficients are 

c v = -3Ar 2/3 ; c v = 4A 2 r 1/3 ; cj = 4 (355) 
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Temperature 



FIG. 14 NbSe 2 phase diagram 
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In the present case, one obtains three different eigenval- 
ues, (de Alameida and Thouless, 1978; Dotsenko, 2001; 
Fischer and Hertz, 1991) 



= 2 (l ± Vl - 4Ar 2 / 3 ) 



(356) 



and — 0. Note that the eigenvalue of Hessian on the 
antisymmetric matrices are degenerate with eigenvalue 
M 1 ) in this case (we will come back later on this eigen- 
value). For A < the solution is unstable due to negative 
/j( 2 ) p or A > 0, both eigenvalues are positive and the so- 
lution is stable. The line A = coincides with the third 
order transition line, hence the non diagonal solution is 
stable when the diagonal is unstable and vise versa. We 
conclude that one of the two RS solutions is stable for any 
value of external parameters (here represented by or and 
r). There still might be a replica asymmetric solution 
with first order transition to it, but this possibility will 
be ruled out within the gaussian approximation and in 
the homogeneous phase without the disorder term 
in the next subsection. Therefore the transition does 
not correspond to RSB. Despite this in the phase with 
nonzero EA order parameter there are Goldstone bosons 
corresponding to in the replica limit of n — > 0. The 
criticality and the zero modes due to disorder (pinning) 
in this phase might lead to great variety of interesting 
phenomena both in statics and dynamics. 

Generalizations and comparison with experi- 
mental irreversibility line 

The glass line resembles typical irreversibility lines in 
both low T c and high T c materials, see Fig. 14, where 
the irreversibility line of NbSe2 is fitted. The theory can 
be generalized to 2D GL model describing thin films or 
very anisotropic layered superconductors. The glass line 



FIG. 15 Melting line and order - disorder lines in layered 
superconductor BSCCO (Beidenkopf et ai, 2007) 



is given in 2D 



2V2 



R-l 



(357) 



An examples of organic superconductor (Shibauchi et al, 
1998) was given in (Li et al, 2006b). The data of BSCCO 
(Beidenkopf et ai, 2005) are compared to the theoretical 
results on Fig. 15. 



4. The disorder distribution moments of the LLL magnetization 

As was discussed in III, the magnetization within LLL 
is proportional to the superfluid density, whose average 
is 



(rP* (r) $ (r)) r = lim - £ (r) * Q (r)) r (358) 

a 

1 ^ 4ttV2 f u k 2 z , 2 _ u , 2 ^ 



U„.a = 2U. 



The variance of the distribution is determined from the 
two thermal averages disorder average: 



(359) 

Within gaussian approximation (Wick contractions) the 
correlators are 



[(^*V>) = 4(u 2 + A 2 ) 



(360) 



46 



Therefore variance of the distribution is given by A. 
This variance determines the width of the magnetization 
loop. In turn, according to the phenomenological Bean 
model (Tinkham, 1996) , the width of the magnetization 
loop is proportional to the critical current. The distri- 
bution of magnetization is not symmetric, as the third 
moment shows: 



8 (u 3 + 3wA 2 + A 3 ) . 



(361) 



It's calculation is more involved. The third irreducible 
cumulant is therefore nonzero: 



(ip*ipY 



- 3 • 2u 



(V>*V>) 2 + 2 (2m) 3 = 8A 3 . (362) 



In analogy to rcf. (Mezard, 1991) one can define "glass 
susceptibility" 



x =<V>»-(V*}<V> = 2(m-A) 



(363) 



useful in description of the "glassy" state. Its variance of 
susceptibility vanishes without RSB: 



X 2 = x 2 - 



(364) 



We return to the replica symmetry breaking after con- 
sidering the crystalline phase. 



C. Gaussian theory of a disordered crystal 

1. Replica symmetric Ansatz in Abrikosov crystal 

In subsection A we used perturbation theory in disor- 
der to assess the basic properties of the vortex crystal. 
However we learned in the previous subsection that cer- 
tain properties like the glass related phenomena cannot 
be captured by perturbation theory and one has to re- 
sort to simplest nonperturbative methods available. In 
the homogeneous phase gaussian approximation in the 
replica symmetric subspace was developed and we now 
generalize to a more complicated crystalline case. This 
is quite analogous to what we did with thermal fluctua- 
tions, so the description contains less details. 
Replica symmetric shift of the free energy 
Within the gaussian approximation the expectation 
values of the fields as well as their propagators serve as 
variational parameters. To implement it, it is convenient 
to shift and "rotate" the fields according to eq.(143): 



ipa (r) = v a ip(r)- 



^2 , Ckexp(ift z )M r )(0k+*^fc)> 



where factor c k = t^t was introduced for convenience 

K l7fcl 

and ipk{x) are the quasimomentum functions. In princi- 
ple the shift as well as fields are replica index dependence. 
However assumption of the unbroken replica symmetry 
means that 



v„ = v, 



(366) 



is the only variational shift parameter. 

To evaluate gaussian energy we first substitute this 
into free energy and write quadratic, cubic and quartic 
parts in fields A and O. The quadratic terms originating 
from the interaction and disorder term are listed below. 
The OO terms coming from the interaction part are: the 
0* a O b term 

[\<P(r)\ 2 [ <Pt(x)c- k c k <p l (r)O a _ k O<l (367) 

Z n Jr Jk,l 

= V ^ J k ^o%o b _ k , 

the 0* a O* b term ^ J k | 7k | O^O b _ k , and finally the O a O b 

term ^- J k | 7k | O k O b _ k . Sum over all four OO terms is 
therefore 

/ ' v 2 (fa + hk\)O a k O b _ k . (368) 
Jk 

Similarly the AA terms sum up to: 

[v 2 ((3 k -\ lk \)AlA b _ k , (369) 

Jk 

while the OA terms cancel. The disorder term con- 
tributes (leading order in n, as usual for replica method): 

V[n / p k J2 A t A -k + I (/3 k -|7k|)5>£A b _ fe ](370) 

2 Jk a J* a,b 

= y I [E(^-|7k|)A^_ fe + ^(/3 k -| 7k |)^_ fe ] 

Z Jk „ , u 



to the A part and 



(365) 



^ 2 [/ £(& + l7k|) 0*0% + £ / (/? k + | 7k |)0£0%] 

1 Jk a a^b Jk 

(371) 

to the O part. The quadratic part of the free energy 
therefore is: 

h = \ I + v " ( 2 & - W) + rv 2 (/3 k - l7k|)]A^ fe 

2Jk a 

+™ 2 E - l7k|) A%A b _ k + [or + v 2 (2/3 k + | 7k |) + rv 2 x 

a^b 

(A + l7k|)]0^0% + ™ 2 E (/?k - |7k|) O a k O b _ k . (372) 

a^b 



47 



There is no linear term and cubic term is not needed since 
its contraction vanishes. The quartic term will be taken 
into account later. We will not need cubic terms within 
the gaussian approximation, while the quartic terms are 
not affected by the shift of fields. We are ready therefore 
to write down the gaussian variational energy. 
Gaussian energy 

Now we describe briefly the contributions to the gaus- 
sian energy. The mean field terms (namely containing 
the shift only with no pairings) are 

f mf = j[nfm- r -n 2 \^ a \\ (373) 

which using eq.(365) takes a form: 



The diagrams are of two kinds. Those including one prop- 
agator and ones which have two propagators from the 
part quartic in fields. The propagators are expectation 
values of pair of fluctuating fields obtained by inverting 
the replica symmetric matrix like in the previous subsec- 
tion. For example for the acoustic mode one gets: 



4n>/2p£ = (A%A a _ k ) = W2; 



4nV2p£ = (A a k A b _ k ) = 



- m 

W2 (Ti£) 2 



(379) 



The integrals of the propagators over k z give: 



fmf = na T v + -Pav 



-n 2 (3 A v 4 



(374) 



The last term can be omitted since the power of n exceeds 
1. Gaussian effective energy in addition to f m f contains 
the Trlog term and the "bubble diagrams". The Trlog 
term comes from free gaussian part, see section III. The 
reference " best gaussian (or quadratic) energy" is defined 
variationally as a quadratic form 



\ J k (e£AtA a _ k +e°O a k O a _ k ) + 
J2[(ri) 2 A a k A b _ k + (jl°) 2 OtO h _ k 



(375) 



and e 



A 
k 



kH2+{^) 2 ;e° = kl/2 + { l £) 2 



(376) 



where ^ , u° , 71^ , Mk > are an variational parameters . We 
assumed no mixing of the A and the O modes, following 
the experience in the clean case and the structure of the 
quadratic part determined in the previous subsection. 

In the following we keep sub leading terms in n since 
they contribute to order n in energy. The Trlog (divided 
by volume) is sum of logarithms of all the eigenvalues. 



1 



trlog 2 3 / 2 7T 2 



1 



2(2tt) 3 

,-: /_4\2 N _i/ 9 (7*k) 



(380) 



((**)- (p£)T 1/a - 

— ~ r / 4ttV2; 
2 (2tt) 3 J kz 

1 r 



■Pk 



Pk 



2tt 



-((/£ 



and similarly for O. 

The contraction of the quadratic parts, after the inte- 
gration and expanded to order n results in: 

k = -\.f trio 6 +nj{p£[»r + « 2 (2/? k - |7k|) - ™ 2 (/3 k - 
bk|)] + rvp^Pk - |7k|)} + n f {p°{a T + v 2 (2f3 k + | 7 k|) 



-rv 2 (A + |7k|)] + rvj%((3 k + | 7 k|)} 



(381) 



For quartic terms coming from two contractions of the 
interaction and the disorder part one obtains, 



fi 



int 



(n-l)log[^-(7i^) 2 ]+log[^ + (n-l) 



/ (p£+p2)fa-i(pf+P?) (382) 

Jk.l 



x (7J^) 2 ] + (n - 1) log[e£ - (/l£) 2 ] + log[ £ ° + („-!) (7l£) 2 ] 
= I f(n - 1)[( M £) 2 - (Ti^) 2 ] 1 / 2 + [(^f + (n - 1) 

+(n - 1)[(^) 2 - (^) 2 ] 1/2 + + (n- 1) 

(377) 

where in the last line integration is over Brillouin zone. 
One observes that the order O (n) terms cancel, while the 
relevant order is: 



+ 



/k,l 

7k 7i 

2/3a 



(rf-rf)(rf-pP) 



and 



his = -\jj{p£+p < Z)^-M+p?) + 



7k7i 





A 



(Pk - Pk) (Pf -P?)\- (p£ + ft) A-, (383) 



(pf+P?) 



7k 7i 

Pa 



Pi 



^-.(Pf-Pf 



tr\ 



n 
2^ 



2[(^) -(^n*+W[(^) -o^r 

+ 2[(^) 2 - (7l£) 2 ^ + (7i°) 2 [(^) 2 - (7l£) 2 ri 

(378) 



jTespectively. Finally we get 

2 

f gauss fmf $2 /tr log .Ant ./dis ' 



(384) 
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2. Solution of the gap equations 
Gap and shift equations 

The gaussian energy is minimized with respect to vari- 
ational parameters. Differentiating with respect to v 2 
one gets the "shift" equation: 

= a T + p A v 2 + J (2/3 k - | 7k |K - r (A - l7k|) 

(385) 

(Pk - Pk) + j[(2& + |7k|)p? - r (Ac + |7k|) (p£ - p£) 

while differentiating with respect to four variational pa- 
rameters in the propagator matrix gap equations are ob- 



tained 








£k = 


a T + 2/3 k w 2 - rv 2 j3 k + (2 - r) 


A k = 


(r 


-i)«V + (i- 


- r) 77k / — 


£ k - 




rw 2 /3 k -r i/pk- 




A k = 


2 f Vl 
Jl PA 




where 








£ k = 


f 

2 


'(/^) 2 +(^) r 




A k = 


f 

2 




^ i 

' Ak= 2 



(387) 



W-(7^) J 



that S„ ~ X™ a T, therefore the coefficients decrease ex- 
ponentially with n. Note that for some integers, for ex- 
ample n = 2, 5,6, /?„ = 0. Retaining only first s modes 
will be called " the s mode approximation" . 

E(k) = E + E lX Pi(k) + ...E nX n p n ( 1 k.)...+ (391) 
£(k) = fy + Eix0iQ*) + -Kx n 0nM- + - 

The expression deviates significantly from the perturba- 
tive one, especially at low temperatures and when the 2D 
case is considered. 

Generalizations and comparison to experi- 
ments 

As we noted already in 2D disorder leads, at least per- 
turbatively, to more profound restructuring of the vortex 
lattice than in 3D. In fact perturbation theory becomes 
invalid. The gaussian methods described above remove 
the difficulty and allow calculation of the order - disor- 
der lime. In this case one does not encounters the "wig- 
gle" but rather a smooth decrease of the order - disorder 
field, when temperature becomes lower. In Fig. 15 the 
ODO line of strongly anisotropic high T c superconductor 
BSCCO is shown. 

One generally observes that there is always off diago- 
nal component in the correlator of the " optical" phonon 
field O. However the off diagonal Edwards - Anderson 
parameter part for a more important low energy excita- 
tion "acoustic" branch appears only below a line quite 
similar to the glass line in the homogeneous phase. 



Solution by the mode expansion 

One can observe that the Ansatz 

(rf) 2 (Mk 4 ) 2 = r/kA; (7l£) 2 - (71k 1 ) 2 = ry k A (388) 

satisfies the gap equations, leading to simpler set for two 
unknown functions and three unknown parameters satis- 
fying eqs. (385), (386), (387) and: 



A 
A 



(1-r) 



(389) 



The equation can be solved by using mode expansion: 

CO 

0k = Y l X n n (ky,0 n {k)= J2 cxp[zk.X] (390) 

"=0 \X\ 2 =na 2 & 

As in IIIC, The integer n determines the distance of a 
points on reciprocal lattice from the origin, and \ = 
exp[— a 2 A /2\ — exp[— 2w/\/3] = 0.0265. One estimates 



D. Replica symmetry breaking 

When thermal fluctuations are significant the efficiency 
of imperfections to pin the vortex matter is generally di- 
minished. This phenomenon is known as "thermal de- 
pinning" . In addition, as we have learned in section III, 
the vortex lattice becomes softer and eventually melts via 
first order transition into the vortex liquid. The inter - 
dependence of pinning, interactions and thermal fluctua- 
tions is very complex and one needs an effective nonper- 
turbative method to evaluate the disorder averages. Such 
a method, using the replica trick was developed initially 
in the theory of spin glasses. It is more difficult to apply 
it in a crystalline phase, so we start from a simpler ho- 
mogeneous phase (the homogeneity might be achieved by 
both the thermal fluctuations and disorder) and return 
to the crystalline phase in the following subsection. 



1. Hierarchical matrices and absence of RSB for the ST C 
disorder in gaussian approximation 

The hierarchical matrices and their Parisi's 
parametrization 
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Experience with very similar models in the theory 
of disordered magnets indicates that solutions of these 
minimization equations are most likely to belong to 
the class of hierarchical matrices, which are comprehen- 
sively described, for example in (Dotsenko, 2001; Fis- 
cher and Hertz, 1991; Mezard, 1991), We limit ourselves 
here to operational knowledge of working with such ma- 
trices contained in Appendix of ref. (Mezard, 1991) 
and collect several rules of using the Parisi's represen- 
tation in Appendix B. General hierarchical matrices u 
are parametrized using the diagonal elements u and the 
Parisi's (monotonically increasing) function u x specifying 
the off diagonal elements with < x < 1 (Parisi, 1980). 
Physically different x represent time scales in the glass 
phase. In particular the Edwards - Anderson (EA) order 
parameter is u x= \ — A > . 

A nonzero value for this order parameter signals that 
the annealed and the quenched averages are different. 
The dynamic properties of such phase are generally quite 
different from those of the non glassy A = phase. In 
particular it is expected to exhibit infinite conductivity 
(Dorsey et al, 1992; Fisher, 1989; Fisher et al, 1991). 
We will refer to this phase as the " ergodic pinned liquid" 
(EPL) distinguished from the "nonergodic pinned liquid" 
(NPL) in which, in addition, the crgodicity is broken. 
Broken ergodicity is related to " replica symmetry break- 
ing" discussed below, however, as we show shortly, in the 
present model of the 5T C disorder and within gaussian 
approximation RSB does not occur. 

In terms of Parisi parameter u and u x the matrix 
equation eq.(324) takes a form: 



Differentiating this equation with respect to x one ob- 
tains; 



+ a T + (4 - 2r) u = 



{u- 2 ) x + 2ru x 



0. 



(392) 



(393) 



Dynamically (see next section), if u x is a constant, pin- 
ning does not results in the multitude of time scales. Cer- 
tain time scale sensitive phenomena like various memory 
effects (Paltiel et al, 2000a,b; Xiao et al, 2002) and the 
responses to "shaking" (Beidenkopf et al, 2005) are ex- 
pected to be different from the case when u x takes mul- 
tiple values. If u x takes a finite different number of n 
values, we call n—1 step RSB. On the other hand, if u x 
is continuous, the continuous replica symmetry breaking 
(RSB) occurs. We show below that within the gaussian 
approximation and the limited disorder model that we 
consider (the 5T C inhomogeneity only) RSB does not oc- 
cur. After that is shown, we can consider the remaining 
problem without using the machinery of hierarchical ma- 
trices. 

Absence of replica symmetry breaking 

In order to show that u x is a constant, it is convenient 
to rewrite the second equation via the matrix a, the ma- 
trix inverse to u: 



x— = 0, 
ax 



(395) 



where we used a set of standard notations in the spin 
glass theory: (Mezard, 1991) 

= M~ (fe) - (fix) = / dxn x ;(396) 

Jo 

[lAx = dy (fj, x - n y ) . 
Jo 

If one is interested in a continuous monotonic part ^ 
0, the only solution of eq.(394) is 



„l/3 



(397) 



{u 2 ) x + 2r(u-% = 0. 



(394) 



Differentiating this again and dropping the nonzero 
derivative again, one further gets a contradiction: 
= . This proves that there are no such monoton- 
ically increasing continuous segments. One can there- 
fore generally have either the replica symmetric solutions, 
namely u x = A or look for a several step - like RSB so- 
lutions (Dotsenko, 2001; Fischer and Hertz, 1991). One 
can show that the constant u x solution is stable. There- 
fore, if a step - like RSB solution exists, it might be only 
an additional local minimum. We explicitly looked for a 
one step solution and found that there is none. 



V. SUMMARY AND PERSPECTIVE 

In this section we summarize and provide references to 
original papers, point out further applications and gen- 
eralizations of results presented here. The bibliography 
of works on the GL theory of the vortex matter is so 
extensive that there, no doubt, many important papers 
and even directions are missed in this brief outline. Some 
of them however can be found in books (Ketterson and 
Song, 1999; Kopnin, 2001; Larkin and Varlamov, 2005; 
Saint- James et al, 1969; Tinkham, 1996) and reviews 
(Blatter et al, 1994; Brandt, 1995; Giamarchi and Bhat- 
tacharya , 2002; Nattermann and Scheidl, 2000). 



A. GL equations. 

Microscopic derivations of the GL equations 

Phenomenological Ginzburg - Landau equations 
(Ginzburg and Landau, 1950) preceded a microscopic 
theory of superconductivity. Soon after the BCS the- 
ory appeared Gorkov and others derived from it the GL 
equations. Derivations and the relations of the GL pa- 
rameters to the microscopic parameters in the BCS the- 
ory are reviewed in the book by Larkin and Varlamov 
(Larkin and Varlamov, 2005) (where extensive bibliogra- 
phy can be found) . The dynamical versions of the theory 
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were derived using several methods and the parameter 7 
in the time dependent GL equation related no the normal 
state conductivity (Larkin and Varlamov, 2005). Most 
of the methods described here can be generalized to the 
case, when the non - dissipative imaginary part of 7 is 
also present. In particular this leads to the Hall current 
(Troy and Dorsey, 1993; Ullah and Dorsey, 1990, 1991) 
and was used to explain the " Hall anomaly" in both low 
T c and high T c superconductors. 

The 5T disorder was introduced phenomenologically 
in statics in (Larkin, 1970). Other coefficients of the GL 
free energy may also have random components (Blatter 
et al, 1994). How these new random variables influence 
the LLL model was discussed in (Li et al, 2006b). 

Anisotropy 

High T c cuprates are layered superconductors which 
can be better described by the Lawrence - Doniach (LD) 
model (Lawrence and Doniach , 1971) than the 3D GL 
model discussed in the present review. The LD model, 
is a version of the GL model with a discretized z coor- 
dinate. However in many cases one can use two sim- 
pler limiting cases. If anisotropy is not very large one 
can use anisotropic 3D GL, eq.(5). The requirement, 
that the GL can be effectively used, therefore limits 
us to optimally doped YBCO7-S and similar materi- 
als for which the anisotropy parameter is not very large: 

7 a = ^Jm*/m* a b = 4 — 8. Effects of layered structure are 

dominant in BSCCO or Tl based compounds (j a > 80) 
and noticeable for cuprates with anisotropy of order 
7a = 50, like LaBaCuO or Hgl223. Anisotropy effec- 
tively reduces dimensionality leading to stronger thermal 
fluctuations according to eq,(19). Very anisotropic lay- 
ered superconductors can be described by 2D GL model 



F = L, 



d 2 r 



, (398) 



which can be approached by the methods presented here. 
For LD model analytical methods become significantly 
more complicated. The gaussian approximation study 
of thermal fluctuations was however performed (Ikcda 
, 1995; Larkin and Varlamov, 2005) and used to ex- 
plained the, so called crossing point of the magnetization 
curves, as well as crossover between the 3D to the 2D be- 
havior(Baraduc et al, 1994; Huh and Finnemore, 2002; 
Junod et al., 1998; Lin and Rosenstein, 2005; Rosenstein 
et al., 2001; Salem-Sugui and Dasilva, 1994; Tesanovic 
et al, 1992) In many simulations this model rather than 
GL was adopted (Ryu et al., 1996; Wilkin and Jensen, 
1997). The GL model can be extended also in direction 
of introducing anisotropy in the a — b plane, like the four 
- fold symmetric anisotropy leading to transition from 
the rhombic lattice to the square lattice (Chang et al, 
1998a, b; Klironomos and Dorsey, 2003; Park and Husc, 
1998; Rosenstein and Knigavko, 1999) observed in many 
high T c and low T c type II superconductors alike (Eskild- 
sen et al., 2001; Li et al, 2006a) . 
Dynamics 



Dynamics of vortex matter can be described by a time 
dependent generalization of the GL equations (Larkin 
and Varlamov, 2005). The bifurcation method presented 
here can be extended to moving vortex lattice (Li et al, 
2004b; Thompson and Hu, 1971). The extension is non- 
trivial since the linear operator appearing in the equa- 
tions is non Hermitian. 

One also can consider thermal transport (Ullah 
and Dorsey, 1990, 1991), for example the Nernst cf- 
fect(Mukerjee and Huse, 2004; Tinh and Rosenstein, 
2009; Ussishkin et al, 2002; Ussishkin, 2003), measured 
recently experiments (Pourret et al, 2006; Wang et al, 
2002, 2006). 



B. Theory of thermal fluctuations in GL model 

Here we briefly list various alternative approaches to 
those described in the present review. It is important 
to mention an unorthodox opinion concerning the very 
nature of the crystalline state and melting transition. 
Although a great variety of recent experiments indicate 
that the transition is first order (for alternative interpre- 
tation see (Nikulov et al, 1995a; Nikulov, 1995b)), some 
authors doubt the existence of a stable solid phase (Kicn- 
appel and Moore , 1997; Moore, 1989, 1992) and therefore 
of the transition all together. 

Functional renormalization group for the LLL 
theory 

While applying the renormalization group (RG) on the 
one loop level, Brezin, Nelson and Thiavillc (Brczin et al., 
1985) found no fixed points of the (functional) RG equa- 
tions and thus concluded that the transition to the solid 
phase, if it exists, is not continuous. The RG method 
therefore cannot provide a quantitative theory of the 
melting transition. It is widely believed however that the 
finite temperature transition exists and is first order (see 
however (Newman and Moore, 1996), who found another 
solution of functional RG equations). 

Large number of components limit 

The GL theory can be generalized (in several differ- 
ent ways) to an N component order parameter field. 
The large N limit of this theory can be computed in 
a way similar to that in the N component scalar mod- 
els widely used in theory of critical phenomena (Itzykson 
and Drouffe, 1991). The most straightforward general- 
ization has been studied in (Affleck and Brczin, 1985) in 
the homogeneous phase leading to a conclusion that there 
is no instability of this state in the 3D GL. However since 
there are other ways to extend the theory to the N com- 
ponent case, it was shown in (Moore et al., 1998) that 
the state in which only one component has a nonzero ex- 
pectation value (similar to the one component Abrikosov 
lattice) is not the ground state of the most straightfor- 
ward generalization. Subsequently it was found (Li et al., 
2004a; Lopatin and Kotliar, 1999) that there exists a gen- 
eralization for which this is in fact the case . 

Diagrams resummation 
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In many body theories one can resum various types 
of diagrams. In fact one can consider Hartree-Fock, 
1/N and even one loop RG as kinds of the diagrams 
rcsummation. Moore and Yeo (Yeo and Moore, 1996a, b, 

2001) and more recently Yeo with his coworkers (Park 
and Yeo, 2008; Yeo et al, 2006) followed a strategy 
used in strongly coupled electron systems to resum all 
the parquet diagrams. The thermal fluctuation in GL 
model had been studied using various analytic methods 
(Koshelev, 1994), but in the vortex liquid region near 
the melting point, or overcooled liquid, non-perturbative 
method must be used, for example, Borel-Pade resum- 
mation method to obtain density density correlation was 
carried out in ref.(Hu et al, 1994). 

Numerical simulations 

The LLL GL model was studied numerically in both 
3D (Sasik and Stroud, 1995) and 2D (Kato and Nagaosa, 
1993; Li and Nattermann , 2003; O'Neill and Moore, 
1993; Tesanovic and Xing, 1991). The melting was found 
to be first order. The results serve as an important check 
on the analytic theory described in this review. In many 
simulations the XY model is employed (Hu et al, 1997; 
Ryu et al, 1996; Wilkin and Jensen, 1997). It is believed 
that results are closely related to that of the Ginzburg - 
Landau model. The methods allows consideration of dis- 
order and dynamics (Chen and Hu, 2003; Nonomura and 
Hu, 2001; Olsson and Teitel, 2001; Olsson, 2007) and fluc- 
tuations of the magnetic field (Nguyen and Sudb0, 1998; 
Sudb0 and Nguyen , 1999). 

Density functional 

The density functional theory is a general method to 
tackle a strongly coupled system. The method crucially 
depends on the choice of the functional. It was applied 
to the GL model by Herbut and Tesanovic (Herbut and 
Tesanovic, 1994) and was employed in (Menon and Das- 
gupta, 1994; Menon et al., 1999; Menon, 2002) to study 
the melting and in (Hu et al, 2005) to the layered sys- 
tems. 

Vortex matter theory 

Elastic moduli were first calculated from the GL model 
by Brandt (Brandt, 1977a,b, 1986) by considering an in- 
finitesimal shift of zeroes of the order parameter. He 
found that the compression and the shear moduli are dis- 
persive. This feature is important in phenomenological 
applications like the Lindemann criterion for both the 
melting and the order - disorder lines (considered dif- 
ferent) (Ertas and Nelson, 1996; Houghton et al., 1989; 
Kierfeld and Vinokur, 2000; Mikitik and Brandt, 2001, 
2003), as well to estimates of the critical current and 
the collective pinning theory (see reviews (Blatter et al, 
1994; Brandt, 1995) and references therein). The disper- 
sion however is ignored in most advanced applications of 
the elasticity theory to statics(Giamarchi and Le Dous- 
sal, 1994, 1995a,b, 1996, 1997; Nattermann and Scheidl, 
2000) or dynamics(Chauve et al, 2000; Giamarchi and 
Le Doussal, 1996, 1998; Giamarchi and Bhattacharya , 

2002) . Recently a phase diagram of strongly type II su- 
perconductors was discussed using a modified elasticity 



theory taking into account dislocations of the vortex lat- 
tice (Dietel and Kleinert, 2006, 2007, 2009) 

Tesanovic and coworkers noted (Tesanovic et al, 1992; 
Tesanovic and Andreev, 1994) a remarkable fact that 
most of the fluctuations effects are just due to condensa- 
tion energy. The lateral correlations part are just around 
2% and therefore can be neglected. The theory explains 
an approximate intersection of the magnetization curves 
and is used to analyze data (Pierson et al, 1995, 1996; 
Pierson and Walls, 1998a; Pierson et al, 1998b; Zhou 
et al, 1993). 

Beyond LLL 

To quantitatively describe vortex matter higher Lan- 
dau levels (HLL) corrections are sometimes required. For 
example in optimally doped YBCO superconductor one 
can establish the LLL scaling for fields above 3T and 
temperature above 70K (see, for example, (Sok et al, 
1995)). A glance at the data however shows that above 
T c the scaling is generally unconvincing: the LLL magne- 
tization is much larger that the experimental one above 
T c . Naively, on the solid side, when the distance from the 
mean field transition line is smaller than the inter - Lan- 
dau level gap, one expects that the higher Landau modes 
can be neglected. More careful examination of the mean 
field solution presented in subsection IIB reveals that a 
weaker condition should be used for a validity test of the 
LLL approximation. However the fluctuation corrections 
involving HLL in strongly fluctuating superconductors 
might be important. Ikeda and collaborators calculated 
the fluctuation spectrum in solid including HLL (Ikeda 
et al, 1990; Ikeda , 1995). In the vortex liquid the HLL 
contribution has been studied by Lawrie (Lawrie , 1994) 
in the framework of the gaussian (Hartree - Fock) approx- 
imation. He found the region of validity of LLL approxi- 
mation. The leading (gaussian) contribution of HLL was 
combined with more refined treatment of the LLL modes 
recently resulting in reasonably good agreement with ex- 
perimental data (Li and Rosenstein, 2003). 

Fluctuations of magnetic field and the dual the- 
ory approach 

Although it was understood that fluctuations of the 
magnetic field in strongly type II superconductors are 
strongly suppressed (Halpcrin et al, 1974; Lobb, 1987), 
they still play an important role when k is not large and 
magnetic field away from H c2 (T) (the situation mostly 
not covered in the present review). The main methods 
are the numerical simulations (Dasgupta and Halpcrin, 
1981; Olsson and Tcitcl, 2003; Sudb0 and Nguyen , 1999) 
and the dual theory approach (Kovner and Rosenstein, 
1992; Kovner et al, 1993; Tesanovic, 1999), which was 
very efficient in describing the Kosterlitz - Thouless tran- 
sition in superconducting thin film and layered materials 
(Ogancsyan et al, 2006). Vortex lines and loops are in- 
terpreted as a signal of "inverted U (1)" or the "magnetic 
flux" symmetry The symmetry is spontaneously broken 
in the normal phase (with photon as a Goldstonc bo- 
son), while restored in the superconductor. Vortices are 
the worldlines of the flux symmetry charges. 
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C. The effects of quenched disorder 

Vortex glass in the frustrated XY model 

The original idea of the vortex glass and the contin- 
uous glass transition exhibiting the glass scaling of con- 
ductivity diverging in the glass phase appeared early in 
the framework of the so called frustrated XY model (the 
gauge glass) (Fisher, 1989; Fisher et ai, 1991; Natter- 
mann and Scheidl, 2000). In this approach one fixes the 
amplitude of the order parameter retaining the magnetic 
field with random component added to the vector poten- 
tial. It was studied by the RG and variational methods 
and has been extensively simulated numerically (Chen 
and Hu, 2003; Chen, 2008; Kawamura, 2003; Nonomura 
and Hu, 2001; Olsson and Tcitcl, 2001; Olsson, 2007).In 
analogy to the theory of spin glass the replica symme- 
try is broken when crossing the GT line. The model ran 
into several problems (see Giamarchi and Bhattacharya 
in ref. (Giamarchi and Bhattacharya , 2002) for a review): 
for finite penetration depth A it has no transition (Bokil 
and Young, 1995) and there was a difficulty to explain 
sharp Bragg peaks observed in experiment at low fields. 

Disordered elastic manifolds. Bragg glass and 
replica symmetry breaking 

To address the last problem another simplified model 
had been proven to be more convenient: the elastic 
medium approach to a collection of interacting line-like 
objects subject to both the pinning potential and the 
thermal bath Langevin force (Cha and Fertig, 1994a, b; 
Dodgson et ai, 2000; Faleski et ai, 1996; Fangohr et ai, 
2001, 2003; Olson et ai, 2001; Reichhardt et ai, 1996, 
2000; van Ottcrlo et ai, 1998). The resulting theory 
was treated again using the gaussian approximation (Gi- 
amarchi and Le Doussal, 1994, 1995a,b, 1996, 1997; Kor- 
shunov, 1993) and RG (Bogner et ai, 2001; Nattermann, 
1990; Nattermann and Scheidl, 2000). The result was 
that in 2 < D < 4 there is a transition to a glassy 
phase in which the replica symmetry is broken follow- 
ing the "hierarchical pattern" (in D = 2 the breaking is 
"one step" ) . The problem of the very fast destruction of 
the vortex lattice by disorder was solved with the vor- 
tex matter being in the replica symmetry broken (RSB) 
phase and it was termed "Bragg glass" (Giamarchi and 
Le Doussal, 1994, 1995a,b, 1996, 1997). A closely re- 
lated approach was developed very recently for both 3D 
and layered superconductors in which effects of disloca- 
tions were incorporated (Dietcl and Kleincrt, 2006, 2007, 
2009). 

Density functional for a disordered systems, 
super symmetry 

Generalized replicated density functional theory 
(Menon, 2002) was also applied resulting in one step 
RSB solution. Although the above approximations to the 
disordered GL theory are very useful in more "fluctu- 
ating" superconductors like BSCCO, a problem arises 
with their application to YBCO at temperature close 
T c (where most of the experiments mentioned above are 
done): vortices are far from being line- like and even their 



cores significantly overlap. As a consequence the behav- 
ior of the dense vortex matter is expected to be different 
from that of a system of line - like vortices and of the 
XY model although the elastic medium approximation 
might still be meaningful (Brandt, 1995). 

One should note the work by Tesanovic and Herbut 
(Tesanovic and Herbut, 1994) on columnar defects in lay- 
ered materials, which utilizes supersymmetry, as an al- 
ternative to replica or dynamics, to incorporate disorder 
non - perturbatively. 

Dynamical approach to disorder in the 
Ginzburg - Landau model 

The statics and the linear response within disordered 
GL model has been discussed in numerous theoretical, 
numerical and experimental works. The glass line was 
first determined, to our knowledge, using the Martin - 
Siggia - Rose dynamical approach in gaussian approxima- 
tion (Dorsey et ai, 1992) and was claimed to be obtained 
using rcsummation of diagram in Kubo formula in (Ikcda 
et ai, 1990). The glass transition line for the AT C disor- 
der was obtained using the replica formalism (within sim- 
ilar gaussian approximation) by Lopatin (Lopatin, 2000) 
and the result is identical to presented in the present re- 
view. He also extended the discussion beyond the gaus- 
sian approximation employing the Cornwall - Jackiw - 
Tomboulis variational method described. This was gen- 
eralized to other types of disorder (the mean free path 
disorder) in ref. (Li et ai, 2006b). The common wisdom 
is that the "replica" symmetry is generally broken in the 
glass (either via " steps" or via " hierarchical" continuous 
process) as in most of the spin glasses theories (Dotsenko, 
2001; Fischer and Hertz, 1991). The divergence of con- 
ductivity on the glass line was obtained in (Roscnstcin 
and Zhuravlev, 2007) (it was assumed in ref. (Dorsey 
et ai, 1992) and linked phenomenologically to the gen- 
eral scaling theory of the vortex glass proposed in ref. 
(Fisher, 1989; Fisher et ai, 1991)). Results are consistent 
with the replica ones presented in the present review) . In 
this work I-V curves and critical current were derived in 
(improved) gaussian approximation and several physical 
questions related to the peak effect addressed. 

Numerical simulation of the disordered 
Ginzburg - Landau model 

The theory can be generalized to the 2D case appro- 
priate for description of thin films or strongly layered su- 
perconductors and compared to experiments. The com- 
parison for organic superconductor k type BEDT—TTF 
(Bel et ai, 2007) and BSC CO (Beidenkopf et ai, 2005) 
of the static glass line is quite satisfactory. There exist, 
to our knowledge, just two Monte Carlo simulations of 
the disordered GL model (Kienappel and Moore , 1997; 
Li and Nattermann , 2003), both in 2D, in which no clear 
irreversibility line was found. However the frustrated XY 
model was recently extensively simulated (Chen and Hu, 
2003; Chen, 2008; Kawamura, 2003; Nonomura and Hu, 
2001; Olsson and Teitel, 2001, 2003; Olsson, 2007) includ- 
ing the glass transition line and I-V curves It shares many 
common features with the GL model although disorder 
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is introduced in a different way, so that it is difficult to 
compare the dependence of pinning. The I-V curves and 
the glass line resemble the Ginzburg - Landau results. 
Finite electric fields 

Finite electric fields (namely transport beyond linear 
response) were also considered analytically in (Blum and 
Moore, 1997) and our result in the clean limit agrees 
with their's. The elastic medium and the vortex dynam- 
ics within the London approximation were discussed be- 
yond linear response in numerous analytic and numer- 
ical works. Although qualitatively the glass lines ob- 
tained here resemble the ones in phenomenological ap- 
proaches based on comparison of disorder strength with 
thermal fluctuations and interaction (Ertas and Nelson, 
1996; Kierfeld and Vinokur, 2000; Mikitik and Brandt, 
2001, 2003; Radzyner et ai, 2002), detailed form is dif- 
ferent. 



D. Other fields of physics 

There are several physical systems in which the meth- 
ods described here, in a slightly modified form, can 
be applied and indeed appeared under different names. 
One area is the superfluidity and the BEC condensate 
physics (ref . (Pethick and Smith, 2008) and references 
therein). Magnetic field is analogous to the rotation of the 
superfluid. One can observe lattice of vortices with prop- 
erties somewhat reminiscent of those of the Abrikosov 
vortices ( Abo-Shaeer et ai, 2001; Baym, 2003; Cooper 
et ai, 2001; Engels et ai, 2002; Madison et ai, 2000; 
Sinova et ai, 2001; Sonin, 2005; Wu et ai, 2007). An- 
other closely relate field is the physics of the 2D electron 
gas in strong magnetic field (Monarkha and Kono, 2004). 
In some cases the problem can be formulated in a way 
similar to the present case with Wigner crystal analogous 
to the Abrikosov liquid (time playing the role of the z di- 
rection of the fluxon), while quenched disorder appears 
in a way similar to the columnar defects in the vortex 
physics. Some aspects of the physics of the liquid crys- 
tal also can be formulated in a form similar to the GL 
equations in magnetic field. 
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VI. APPENDICES 

A. Integrals of products of the quasimomentum 
eigenfunctions 

In this appendix a method to calculate space averages 
of products of the quasi - momentum eigenfunctions in 
both static and the dynamic cases. 



1. Rhombic lattice quasimomentum functions 

Let us specialize to a rhombic lattice with following 
bases of the direct and the reciprocal lattices (see Fig. 3 
for definition of the angle 9 and the lattice spacing a$, 
subject to the flux quantization relation, eq.(63)): 

di = d e (l,0), d 2 = d fl (i,itantf); (399) 

di = d e (itan#,-i), d 2 = (0,dj). 

We use here the "LLL" unit of magnetic length. 

We start with static LLL functions for an arbitrary 
rhombic lattice: 



dg 



E 



x e 



-Hv-k*-%i) 2 



(400) 



To include higher LL corrections, it is conve- 
nient to use rising and lowering operators intro- 
duced in eq.(90) to work with the HLL func- 
tions, ot = (2)~ 1/2 (-id x +d y -y) , a = 
— (2) ' (—id x + d y + y). The following formula will be 
frequently used. If ip is an LLL function, then 



^*a +N f (x, y) = 2-^{-id x + d y ) n V *f (x, y) . (401) 

2. The basic Fourier transform formulas 
Product of two functions 
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It can be verified by direct calculation of gaussian in- and the product takes a form: 
tegrals that 



^ 2 x 



J <Pk( r )<Pi( r )<Pv( r )<Pv( r ) = (2tt) 2 
L{r)^{vy^=A^ 5(JC-k-K)F(k,K); V J [k- 1- (k' - 1') - A] / [k, 1, k', 1', A] (411) 

F(k,K)=e-^ +i ^- E ^ L+k ^\ (402) / [k, 1, k', 1', A] = ]T e^x 

with decomposition of arbitrary momentum JC into its K,K _K + A 

"rational part" k, which belongs to the Brillouin zone and e l If K l~% K 'i 2 +( k ^- k L) lc v- lc ^( l y- l 'y)+ l v( k - l ) x - l 'y{ k ' 
an "integer part" K, belonging to the reciprocal lattice 

where / [k, 1, k', 1', A] = if k - 1- (k' - 1') - A ^ 0. 
/C = k + K, K =K 1 d 1 +K 2 d 2 . (403) The last exponent in function / [k, 1, k', 1', A] can be also 

Its inverse Fourier transform, rearranged as 

K 2 7T? 

• h r -^--A 1 (2X 1 + A 1 )+z(k-k') (412) 

</>(r)y>*(r + k) =e lxk *J2e- lK - r F(k,K), (404) 2 2 y ' 

K AJC+iJC x A y + i (l y - l' y ) (k - l) x + il' y A x . 

where ki — kj , can be generalized into Using 

^(r+lV(r + k) = e < ( I+, »N fc - , >. x (405) e jXq = cell ^ <5 (q - K) , 

X=nidi+n 2 d2 K 



^ e -iK-(r+i) + ^(A-?+/{-i)-^^^+i(fe-0 x /C s - — 
K 



where ceZZ is the volume of the unit cell and in our units 
is equal to 2ir, one obtains the Poisson resummation re- 
with JC =K + k 1. This in turn provides a very useful lation 
product representation: 

riM*.(r) = £ e----^ (406) £ , (K) = J_ | ^ cxp (<x . q) , (q) > 



Xe'f K l~' E^ES +ik x )Cy—i)Ca;ly+ily(k—l) a . 

The four - point vertex function 



K "X 

Using Poisson resummation, one rewrites the sum as 



±ne jour - point verzex junczion , i(x+?x(k-k')) 2 -*x.(k-i) 

The relation above used twice gives the following ex- ■' I ' ' ' ' J 
pression for the four point vertex:in the quasi - momen- x 
turn representation: e ^(ky-k' y )(k x -l x )+ly^x-^ A ?l _ (413) 

/ lpI.(y)lpi(y)lpI (r)^k' (r) = (27r) 2 £ <5 [/C - JC'] (407) 3. Calculation of the /3 k , 7 k functions and their small 

vT 1-- mnmontnm ovn^ncinn 



xe -T-+*[f^i 2 -f^; 2 +(^- fc x)'c H -'^(i ! ,-^)+i H ('=-o :c -';(fe'-'') 



AC momentum expansion 

One often encounters the following space averages: 



The delta function S [JC — JC'] implies that 

^^MVA-^^^) 2 ^). (414) 

K 1 + ki-h = K' 1 + h{-l' 1 . (408) 

/3k = /?k an d 7k = 7jf- Using formulas of the previous 
As < kiji,k[j[ < 1, we have only three possible inte- subsection, one can write 
ger values for each coordinate: 

¥>Vk = N/2 r— (~id x + d y f = 
fci-Ji-fci+Ji = A! = 0,1,-1, (409) 1 V ' 

fc 2 -/ 2 -4 + ^ = A 2 = 0,l,-1, ^^Y( Z * + z Q ) N e^>*F*(KQy, (415) 

which require /Ti^ — K[ 2 = 0, —1, 1. Thus ^ 

_ ' wf = (<*>Vk )* = 



2 w / 2 v / M 

iV 



[k — 1— (k' — 1') — A] ; Q 
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Therefore 



N 



K= y (^x) 



iV e _2|l_ ik .x 



and 



fc = E 



g-^-ik-X 



Similarly 7^ can be obtained, for 7 k , we have 



7k = e 



E< 



-5 4Z k ZX 



(416) 



(417) 



(418) 



When the symmetry is higher, the expressions simplify. 
Using the six - fold symmetry of the hexagonal lattice, 

X' x + iX' y = e w (X x + iX y ) ,6=^1; (424) 

the sum eq.(419) transforms into 

S (N, M) — S (N, M) e ln0 . (425) 

Thus S (N, M) = if N ^ 6j. Using S (2, 0) = S (4, 0) = 
S (2,2) = 0, one obtains several relations of different 
sums to simplify expansion of /3k 



a V- -* f1 k2X2 fc4 ^\ 



64 



(426) 



The above formula is valid for any lattice structure. 
Small momentum expansion of the /3k, Jk func- 
tion for the general rhombic lattice 

Consider the sum 



5(7V,M) = ^e-^^|X| M 



for any integers N, M. Due to reflection symmetry 

^e-^Jf„ b = (420) 
x 

for l\, I2 integers, and h + h odd integer. For small k: 

(ik-X)' 



a-E^+E^) 



(421) 



- /3A-E 



(^ll + kyXy) 



E< 



24 



+ 0(k 6 



Similarly for the function 7^ can be expanded for small 
k 2 



7fe =e **-*» ^E e ^t 1 + E 



Q 



(422) 



so that 



00 *2,l 2/ 

h*l = e-E^(l + E^ 1/2 (423) 

Q 1=1 ^ 

00 Jll' *2l' 



Q' 



Small momentum expansion of the /3k, Jk func- 
tion for hexagonal lattice 



Similarly 



M^All-y + yRO^ 6 ), (427) 



(419) and its phase 6 k has an expansion 



7k 
l7k| 



= 1 - ik x k y + O (fc 4 ) ; 6 k = -k x k y + O (fc 4 ) . (428) 



In terms of z* it is an analytic function: 



7fe = e 



/.A.i: h.- n 



5eZ/ duality relation 

If the lattice is self-dual, one can prove 



E*v*=/? A . 



(430) 



Thus 



/?k = /?A-^ 2 + gE^ e "^- ( 431 ) 

X 

Using the expansion for /3k, 7k, one can obtain the su- 
persoft acoustic phone spectrum: 



- /9a + 2/3 k - | 7 k| = E ~ ^) fc4 + ( fc6 ) • 

(432) 

T/ie small momentum expansion of the vertex 
function 

For momentum, k, 1, k' are not too big to have Umk- 
lapp process 

J^(r)v>i(r)p5(r)^(r) = (2^) 2 <5 [k - 1-k'] x 

+*( fe -- fc xK«-*Mi ! ,-i;)+*(i H )(fc-z) x , (433) 

K 
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where /C=K + k — 1. In this case, we define 
/ [k, 1, k', 1', A] = [k, l'|l, k'] and which has small 
momentum expansion: 



[k,0|l,k'] =/3 A (l 



P + k' 2 



+ ^ {kxky l x l y k' x ky)) 
(434) 

Another useful identity 

Any sixfold (Dq) symmetric function F(k) (namely a 
function satisfying -F(k) = F(k ), where k, k is related 
by a ^ rotation) obeys: 



L 



^(k)7k7k,i = / F(k) |7k| 2 ,7k.i =< <pt<p*_ k <p-m > 

PA J k 



(435) 

This identity can be proved by using the fact 
(f k -F(k)7k7k,i) / 7i is a analytic function of z\ and is a 
periodic function of reciprocal lattice vectors, i.g. 



1 — » 1+midi + m 2 d 2 



(436) 



the function is unchanged. The only solution for a ana- 
lytic function with this property is a constant. 



/ k F (k)7k7k,i 



const. 
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(437) 



The constant can be determined by setting 1 = 0. 



B. Parisi algebra for hierarchial matrices 

In this appendix we collect without derivation the for- 
mulas used in calculation of disorder properties. Deriva- 
tion can be found in ref. (Mezard, 1991). Inverse matrix 
has the following Parisi parameters: 



m 



to. 



= . 1 a-/' 

TO— < TO > J Q 



du 



TO 



to 



TO — < TO > ' 
1 



u m— < to > — [m\ u 

(438) 

[ml. 



+ 



to — < to > x [to— < to > — [m\ x \ 
dv [toL toq 



I 

Jo 



i Q it to — < TO > — |roj„ to— < TO > 
Square of matrix can be treated similarly: 



) 



a.b 



(to) 2 - < (to-j;) 2 >; 



(439) 



(to 2 )^ = 2 (to— < to >) m x — / 

Jo 



dv(m x — m v ) 2 
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